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Abstract 

We  present  a  new  stability  analysis  for  hybrid  legged 
locomotion  systems  based  on  the  “symmetric”  fac¬ 
torization  of  return  maps.  We  apply  this  analysis 
to  2  and  3  degree  of  freedom  (DOF)  models  of 
the  Spring  Loaded  Inverted  Pendulum  (SLIP)  with 
different  leg  recirculation  strategies.  Despite  the 
non-integrability  of  the  SLIP  dynamics,  we  obtain  a 
necessary  condition  for  asymptotic  stability  (and  a 
sufficient  condition  for  instability)  at  a  fixed  point, 
formulated  as  an  exact  algebraic  expression  in  the 
physical  parameters.  We  use  this  expression  to  study 
a  variety  of  2  DOF  SLIP  models  that  have  previously 
been  posited  as  low  dimensional  representations  of 
running,  focusing  on  the  sensory  “cost”  required  to 
achieve  “fast”  transients  as  measured  by  the  degree  of 
singularity  of  the  linearized  dynamics.  We  introduce 
a  new  3  DOF  SLIP  model  with  pitching  dynamics 
whose  stability  properties,  revealed  by  this  analysis, 
provide  for  the  first  time  the  beginnings  of  a  formal 
explanation  for  the  surprisingly  stable  gaits  of  the 
open  loop  controlled  robot,  RHex. 

Keywords  -  legged  locomotion,  hybrid  system, 
return  map,  Spring  Loaded  Inverted  Pendulum, 
stability,  time-reversal  symmetry 

1  Introduction 

This  paper  introduces  a  new  formalism  for  studying 
the  stability  of  legged  locomotion  gaits  and  other  pe¬ 
riodic  dynamically  dexterous  robotic  tasks.  We  are 
motivated  in  part  by  the  need  to  explain  and  control 
the  remarkable  performance  of  RHex,  an  autonomous 


hexapedal  running  machine  whose  introduction  broke 
all  prior  published  records  for  speed,  specific  resis¬ 
tance,  and  mobility  over  broken  terrain  [1],  Powered 
by  only  six  actuators,  located  at  the  “hips”  to  drive 
each  of  its  six  passively  compliant  legs,  in  the  manner 
of  single-spoked  rimless  wheels,  RHex’s  locomotion  is 
excited  by  a  single  periodic  “clock”  signal  split  into 
phase  and  anti-phase  copies  for  coordinating  its  alter¬ 
nating  tripod  gait.  A  simple  PD  controller  at  each  hip 
motor  in  a  given  tripod  forces  its  leg  to  track  the  alter¬ 
nately  fast  and  slow  clock  reference  signal  correspond¬ 
ing  to  presumed  stance  and  swing  phases.  Experi¬ 
mentally,  RHex’s  performance  at  various  speeds  over 
various  terrains  is  strongly  dependent  upon  the  par¬ 
ticular  values  of  the  clock  parameters,  and,  as  is  typi¬ 
cal  within  the  feedforward  control  paradigm,  each  new 
situation  demands  its  own  carefully  tuned  parameter 
set.  Better  analytical  understanding  of  the  relation¬ 
ship  between  clock  signal  and  steady  state  gait  should 
dramatically  simplify  the  frequently  lengthy  empiri¬ 
cal  parameter  tuning  exercises  presently  required  to 
achieve  high  performance  gaits. 

1.1  The  SLIP  Model  as  a  Template  for 
RHex 

A  complete  account  of  the  relationship  between 
RHex’s  clock  signal  and  steady  state  gait  in  even  the 
simplest  case  would  entail  insight  into  the  steady  state 
properties  of  an  under-actuated  high  degree  of  freedom 
(DOF)  hybrid  mechanical  system  whose  Lagrangian 
dynamics  switches  among  a  set  of  26  possible  holo- 
nomically  constrained  models  depending  upon  which 
feet  are  in  contact  with  the  ground.  Fortunately,  a 
growing  body  of  simulation  study  and  empirical  evi¬ 
dence  [2]  suggests  that  RHex,  when  properly  tuned, 
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exhibits  sagittal  plane  stance  behavior  well  approxi¬ 
mated  by  the  two  degree  of  freedom  SLIP.  This  comes 
in  addition  to  the  well-established  fact  that  the  cen¬ 
ters  of  mass  of  running  animals  [3]  and  humans  [4]  ap¬ 
proximately  follow  the  dynamics  of  the  sagittal  plane 
SLIP.  A  general  framework  for  “anchoring”  the  SLIP- 
“template”  mechanics  in  the  far  more  elaborate  mor¬ 
phologies  of  real  animals’  bodies  has  been  introduced 
in  [5].  Briefly,  given  a  high-dimensional  dynamical 
system  -  the  “anchor”  -  that  is  believed  to  be  a  rea¬ 
sonably  accurate  model  of  an  animal  or  robot,  a  “tem¬ 
plate”  is  a  low-dimensional  dynamical  system  whose 
steady  state  encodes  the  task  and  is  conjugate  to  the 
restriction  dynamics  of  the  anchor  on  an  attracting 
invariant  submanifold. 

In  general,  both  the  anchoring  as  well  as  the  con¬ 
trol  of  the  SLIP  template  seem  to  demand  sensing, 
actuation,  and  computation  that  may  be  unrealistic 
relative  to  the  resources  that  animals  and  practical 
robots  might  possess.  Indeed,  a  hierarchical  controller 
[6]  for  a  RHex-like  simulation  model  programmed  in 
SimSect  [7]  that  enforces  both  the  anchoring  as  well  as 
the  template  control  relies  on  sophisticated  full  state 
feedback.  A  part  of  the  sensor-suite  necessary  to  im¬ 
plement  this  feedback  control  has  only  recently  been 
installed  on  the  robot  [8]  and  it  is  currently  unknown 
whether  the  stabilizing  effect  of  this  controller  seen  in 
simulation  will  persist  in  the  presence  of  unavoidable 
sensor  noise.  This  motivates  the  question:  is  it  pos¬ 
sible  to  implement  the  template-anchor  paradigm  [5] 
with  sensor-cheap,  low-bandwidth  robotic  controllers? 
In  this  paper  we  address  that  part  of  the  above  ques¬ 
tion  concerned  with  template  control.  Namely  given 
that  a  SLIP-anchoring  mechanism  is  present,  either  by 
deliberate  design  or  by  the  interaction  of  the  controlled 
robot  with  its  environment,  can  the  stability  and  per¬ 
formance  of  the  controlled  template  be  assessed  me¬ 
thodically  (beyond  empirical  or  numerical  study),  for 
example,  as  a  function  of  the  cost  of  the  sensory  feed¬ 
back  required? 

1.2  Output  feedback  stabilization  in 
the  SLIP  model 

The  SLIP  model  is  a  hybrid  dynamical  system  formed 
by  the  composition  of  leg-body  stance  dynamics  with 
ballistic  body  flight  dynamics.  Control  takes  place 
during  the  flight  phase,  where  the  leg  angle  is  set  for 
the  next  touchdown  event.  The  two  degree  of  free¬ 
dom  SLIP  model  provides  a  ubiquitous  description  of 
biological  runners  in  the  sagittal  plane  [3]  and,  as  men¬ 
tioned  above,  a  broadly  useful  prescription  for  legged 
robot  runners  such  as  RHex  [9,  1,  2]  as  well.  The 


closely  related  three  degree  of  freedom  Lateral  Leg 
Spring  (LLS)  model,  has  been  recently  identified  as  a 
candidate  template  for  cockroach  running  in  the  hor¬ 
izontal  plane  [10,  11]  and  seems  likely  to  be  relevant 
for  RHex  as  well  [1]. 

However,  the  limitations  of  the  two  degree  of  freedom 
SLIP  model  (no  pitching  dynamics,  no  lateral  dynam¬ 
ics)  and  the  three  degree  of  freedom  LLS  model  (fail¬ 
ure  to  reproduce  some  aspects  of  animal  data  [12]) 
show  that  far  more  sophisticated  models  will  be  re¬ 
quired  to  capture  more  salient  features  of  the  anchor. 
In  particular,  a  literal  template  of  RHex,  i.  e.  a  model 
whose  dynamics  represents  the  restriction  dynamics 
of  an  attracting  invariant  submanifold  in  RHex,  must 
include  a  source  of  dissipation  as  well  as  hip  torques. 
Despite  these  shortcomings,  the  two  degree  of  freedom 
SLIP  and  its  extension  to  three  degrees  of  freedom  (in¬ 
troduction  of  pitch  dynamics)  are  sufficiently  well  mo¬ 
tivated  by  prior  literature,  sufficiently  mathematically 
challenging  (due  to  their  non-integrable  nature)  and 
their  analysis  sufficiently  revealing  of  RHex-like  prop¬ 
erties  (as  is  shown  below)  as  to  motivate  our  exclusive 
focus  on  them  in  this  paper. 

The  stability  properties  of  these  hybrid  dynamical  sys¬ 
tems  can  be  assessed  by  a  Poincare  or  return  map  R 
acting  on  a  (reduced)  Poincare  section  X : 

R  :  X  X  .  (1) 

In  legged  locomotion,  the  iterates  of  this  return  map 
R  -  the  function  relating  the  body  state  at  a  periodi¬ 
cally  (at  each  stride)  occurring  event  -  summarizes  all 
properties  relevant  to  the  goal  of  translating  the  body 
center  of  mass.  The  return  map  arises  in  general  from 
a  controlled  plant  model 

x(k  + 1)  =  A(x(k),u(k)) 

y{k)  =  C(x(k))  (2) 

where  the  discrete  time  control  input  variable,  u(k), 
represents  the  consequences  at  the  integrated  stride- 
by-stride  level  of  controlled  influences  imposed  over 
continuous  time  within  stance  or  flight.  In  this  pa¬ 
per,  physically  motivated  assumptions  (listed  in  Sec¬ 
tion  2.4.1)  that  we  impose  upon  the  allowable  contin¬ 
uous  time  influences  turn  out  to  yield  a  discrete  time 
representative,  u,  that  implicitly  determines  the  flight 
time  for  the  ballistic  phase  of  the  body  at  each  stride. 
When  the  continuous  time  physical  influences  imposed 
within  a  given  stride  are  determined  according  to  state 
information  gathered  from  the  available  observations 
of  the  previous  stride,  we  have  effectively  introduced 
a  discrete  time  feedback  policy, 

u(k)  =  H{y{k))  (3) 
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whose  closed  loop  yields  (1),  R(x )  =  A(x,H  o  C(x)). 
The  controlled  plant  model  for  SLIP  systems  is  spec¬ 
ified  in  Section  2.4.3. 

In  this  paper  we  will  confine  our  study  exclusively  to 
such  time  invariant  output  feedback  laws,  H  (3)  for 
two  allied  reasons.  First,  this  restriction  focuses  at¬ 
tention  on  the  key  role  played  by  the  output  function, 
C  (2),  variations  of  which  we  will  use  to  model  sensor 
limitations  of  the  underlying  physical  system  repre¬ 
sented  by  the  SLIP  model.  Second,  since  u  models 
the  influence  of  flight  phase  duration  (implicitly  by 
specifying  the  leg  angle  trajectory),  this  restriction  to 
time  invariant  output  feedback,  H  (3),  models  the  leg 
recirculation  policies  that  have  so  rightly  captured  the 
attention  of  the  legged  locomotion  community  in  re¬ 
cent  years. 

The  surprising  discovery  of  “self-stable”  legged  loco¬ 
motion  -  first  in  the  closely  related  LLS  model  [11]; 
subsequently  in  the  SLIP  itself  [13,  14]  -  demands  a 
more  systematic  account  of  what  is  meant  by  the  term 
“self.”  In  these  studies,  the  duration  of  flight  phase  is 
determined  by  a  fixed  leg  angle  policy,  and  “self’  con¬ 
notes  the  apparent  absence  of  active  sensors.  Recently, 
a  more  elaborate  state-dependent  leg  retraction  policy 
has  been  shown  numerically  to  inherit  the  stability 
properties  of  the  fixed  touchdown  angle  policy  while 
increasing  the  basin  of  the  stable  gait  [15].  On  the 
other  hand,  a  recirculation  policy  that  initiates  after 
leg  liftoff  a  constant  angular  velocity  until  leg  touch¬ 
down  can  induce  neutral  stability  [16].  These  appar¬ 
ently  slightly  varied  policies  mask  significant  variation 
in  cost  and  effort  depending  upon  how  the  sensor  suite 
might  be  implemented  in  practice.  We  seek  to  shed 
greater  light  on  when  a  more  or  less  clever  leg  recircu¬ 
lation  strategy  can  make  a  difference  in  the  quality  of 
gait  stability  (e.g.,  faster  transients,  larger  basin)  as  a 
function  of  the  “cost”  of  sensory  data. 

Of  course,  real  sensors  are  not  implemented  in  these 
templates  at  all  but  in  physical  machines.  Empirically, 
it  is  abundantly  clear  that  the  leg  swing  policy  plays 
a  central  role  in  the  gait  quality  of  physically  useful 
machines  such  as  RHex  [1].  Leg  recirculation  strate¬ 
gies  have  been  shown  numerically  to  play  a  key  role 
in  the  gait  quality  of  independent  locomotion  models 
inspired  by  quadrupedal  animal  trotters  [17]  and  gal¬ 
lopers  [18].  When  the  SLIP  template  is  anchored  ac¬ 
tively  [6]  then  its  stability  properties  determine  those 
of  the  anchor  by  definition,  hence  insight  into  how  to 
tune  the  quality  of  SLIP  gaits  transfers  directly  over  to 
the  physical  machine  of  interest.  However,  in  the  ab¬ 
sence  of  complete  state  feedback,  the  correspondence 
between  the  template  and  the  behavior  of  a  complex 
system  that  shows  empirical  evidence  of  anchoring  it 


is  not  at  all  clear. 

Lacking  formal  results  bearing  on  this  issue,  we  find  it 
useful  to  introduce  terminology  summarizing  the  fol¬ 
lowing  intuitive  distinction.  We  will  say  that  the  cor¬ 
respondence  is  “descriptive”  if  properties  observed  in 
the  complex  model  are  also  observed  in  the  template 
model  fitted  to  it.  We  will  say  that  the  correspon¬ 
dence  is  “prescriptive”  if  design  parameter  settings  in 
the  complex  model  indeed  anchor  and  yield  as  well  the 
same  properties  they  produce  in  the  template  model. 

1.3  Contribution  of  this  paper 

Notwithstanding  its  apparent  simplicity,  the  SLIP 
model  is  non-integrable:  the  stance  phase  trajectory 
cannot  be  written  down  in  closed  form  [19].  This  has 
motivated  authors  who  seek  insight  more  systematic 
than  numerical  simulation  can  provide  to  develop  var¬ 
ious  physically  motivated  closed  form  approximations 
to  R  instead  [20,  21,  22].  In  contrast,  here  we  observe 
that  while  R  cannot  be  written  in  closed  form,  cer¬ 
tain  physically  reasonable  assumptions  (listed  in  Sec¬ 
tion  2.4.1,  below)  imply  that  the  determinant  of  its 
Jacobian  at  a  symmetric  fixed  point  (to  be  defined  in 
Section  2.3)  of  R  can  be  so  expressed.  The  central  con¬ 
tributions  of  this  paper  arising  from  that  observation 
include: 

1.  A  new  analytical  framework  based  on  a  “symmet¬ 
ric”  factorization  of  the  return  map  R ,  in  terms  of 
its  non-hybrid  components  that  yields  the  closed 
form  expression  of  the  determinant  at  a  symmet¬ 
ric  fixed  point  of  R  (Section  3).  Necessary  condi¬ 
tions  for  asymptotic  stability,  sufficient  conditions 
for  instability,  and  conditions  equivalent  to  neu¬ 
tral  stability  of  the  closed  loop  map,  i?,  follow. 

2.  Closed  form  conditions  on  HoC  yielding  rigorous 
statements  concerning  the  sensory  “cost”  of  con¬ 
trol  in  both  the  2  DOF  and  3  DOF  settings  that 
cannot  be  established  by  mere  numerical  study, 
as  follows: 

(a)  2  DOF  SLIP  models:  any  control  with  fast 
transients  (“singular”  control:  the  Jacobian 
of  the  closed  loop  return  map  is  globally  sin¬ 
gular)  requires  velocity  sensing  and  is  there¬ 
fore  “costly”  (Section  3.3.1). 

(b)  3  DOF  SLIP  models:  SLIP  models  that  have 
only  non-inertial  (body  frame)  sensors  avail¬ 
able  cannot  implement  singular  control  (Sec¬ 
tion  3.4.1). 

3.  A  new  3  DOF  SLIP  model  based  upon  the  RHex 
gait  generator  [1]  which,  when  subject  to  the 
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factored  stability  analysis,  imposes  for  the  first 
time  analytical  conditions  on  RHex’s  “feedfor¬ 
ward  clock”  parameters  [1]  necessary  for  the  sta¬ 
bility  (and  sufficient  for  instability)  of  the  result¬ 
ing  3  DOF  gait  (Section  3.4.2). 

4.  A  preliminary  numerical  study  in  Section  4  sug¬ 
gesting  that  our  3  DOF  SLIP  model  is  a  good 
descriptive  model  of  a  RHex-like  24  DOF  model 
(programmed  in  SimSect  [7]). 

In  Section  2  we  preface  this  analysis  by  introducing 
the  terminology  and  notation  for  hybrid  systems  to 
be  used  subsequently,  followed  by  a  review  of  how  re¬ 
versibility  symmetries  can  replace  the  symplectic  sym¬ 
metry  in  Liouville’s  theorem  (see  e.g.  [23]),  which  does 
not  generally  apply  to  hybrid  systems.  We  then  de¬ 
velop  the  consequences  of  these  observations  in  Section 
3  along  the  lines  described  above,  present  a  prelimi¬ 
nary  numerical  study  in  Section  4,  and  close  with  some 
brief  concluding  remarks  in  Section  5. 

2  Theoretical  framework  and 
modeling  assumptions 

In  section  2.1  we  introduce  the  terminology  of  hybrid 
dynamical  systems  and  provide  some  intuition  con¬ 
cerning  the  machinery  used  to  trim  away  the  awk¬ 
ward  and  inessential  details  of  our  hybrid  model  to 
yield  a  conventional  discrete  dynamical  control  sys¬ 
tem  (2)  whose  closed  loop  properties  (1)  represent  the 
formal  object  of  study.  Having  established  a  notation 
for  (hybrid)  dynamical  systems,  Liouville’s  theorem,  a 
key  tool  in  the  present  study,  can  be  stated  formally 
in  the  next  section,  2.2.  Then  an  analogue  of  the  local 
form  of  Liouville’s  theorem  for  discrete  maps  derived 
from  hybrid  systems  will  be  established  in  Section  2.3. 
Section  2.4  formally  defines  the  SLIP  system  with  its 
hybrid  components  as  well  as  its  Poincare  section  and 
discrete  time  return  map. 

2.1  Preliminary  definitions  and  mod¬ 
eling  considerations  of  hybrid  dy¬ 
namical  systems 

Models  of  legged  locomotion  are  characterized  by  dis¬ 
tinct  phases,  notably,  stance  and  flight.  Formally,  the 
dynamics  cannot  be  described  by  a  single  flow,  but 
require  a  collection  of  continuous  flows  and  discrete 
transformations  governing  their  transitions.  The  re¬ 
sulting  model  is  called  a  “hybrid”  system.  This  section 
makes  the  notion  of  a  hybrid  system  more  precise  by 
adapting  the  definitions  in  [24]  to  the  present  setting. 


Let  X  be  a  finite  index  set  and  Xa ,  a,  £  X  with 
dim(Aa)  =  2 N  a  collection  of  open  Euclidean  domains 
(charts).  Assume  a  mechanical  system  whose  time 
evolution  is  described  by  holonomically  constrained 
autonomous  conservative  vector  fields  fa,  with  con¬ 
figuration  space  variables  q:  x  =  fa{x)  with  x  = 
( q  q)T  £  Xa.  Transitions  from  one  vector  field  fa 
to  another  vector  field  fp  are  governed  by  thresh¬ 
old  functions  hf^  which  specify  an  event  at  their  zero¬ 
crossing.  The  threshold  functions  h ^  can  depend  on 
the  initial  condition  xq  =  x(t  =  0)  €  Xa,  time  t,  and 
the  current  state  ff(x o)  =:  x{t)}  We  restrict  our¬ 
selves  to  hybrid  systems  where  for  each  chart  there 
is  only  one  threshold  function  /),((;  hence  the  upper 
index  (3  will  be  dropped  from  now  on.  We  also  re¬ 
set  the  time  to  zero  at  each  chart  transition.  The 
end  time  of  the  evolution  on  chart  Xa  is  uniquely  de¬ 
fined  by  ta( x0)  =  mint>0{f  :  ha(fta(x0),x0R)  =  0}. 
The  equation  ha(f^(xo),Xo,t)  =  0  will  be  referred  to 
as  the  threshold  equation.  Switching  between  charts 
is  effected  by  transition  mappings  T@  with  domains 
in  Xa  and  ranges  in  Xp.  The  flow  map  Fa  for  the 
ath  vectorfield  is  defined  via  the  implicit  function,  ta, 
Fa:x 0w  /a“(xo)(a: 0).2 

In  this  paper,  as  in  many  settings  of  hybrid  dynamical 
systems,  we  are  interested  in  the  attractive  behavior 
of  distinguished  orbits  whose  appropriate  projections 
are  periodic.  By  “periodic”  we  mean  that  the  dis¬ 
tinguished  orbit  is  defined  on  a  recurring  sequence  of 
charts  along  which  the  projected  flow  yields  a  return 
to  the  same  projected  initial  condition.  An  “appropri¬ 
ate”  projection  strips  away  variables  whose  values  are 
not  descriptive  of  the  locomotion  task  -  here,  the  con¬ 
served  total  mechanical  energy  along  with  the  cyclic 
variable  of  elapsed  distance.  Similarly,  “attractive  be¬ 
havior”  denotes  the  asymptotic  properties  of  projected 
orbits  relative  to  the  projection  of  the  distinguished  or¬ 
bit.  These  slight  variants  of  the  traditional  Poincare 
analysis  of  conventional  dynamical  systems  theory  will 
all  be  introduced  formally  in  the  next  section,  and  will 
be  seen  to  yield  a  stride  map 

S  =  S2oS1  (4) 

whose  projection  (along  with  those  of  its  factors,  Sa) 
that  we  will  denote  R  (along  with  the  corresponding 
factors,  Ra )  captures  as  a  discrete  time  iterated  dy¬ 
namical  system  the  locomotion  relevant  behavior  of 

1  Note  that  this  is  more  general  than  the  definition  in  [24], 
where  hf  only  depends  on  //,  (xq  ) .  This  added  generality  is 
required  because  we  wish  to  study  the  effects  of  leg  recirculation 
strategies  -  reference  trajectories  parametrized  by  time  that  are 
triggered  by  a  liftoff  transition  -  see  e.g.  equation  (47). 

2Note  that  Fa  is  not  the  usual  constant-time  flow  map  of 
dynamical  systems  theory  o);  rather  the  time  varies  de¬ 
pending  upon  the  initial  data  xq. 


4 


our  hybrid  dynamical  system  analogous  to  a  Poincare 
map. 

2.2  Liouville’s  theorem  and  stability 

Informally,  Liouville’s  theorem  states  that  volume  in 
phase  space  of  a  holonomically  constrained  conserva¬ 
tive  dynamical  system  described  by  a  single  Hamilto¬ 
nian  flow  is  preserved,  i.e.  a  set  of  initial  conditions 
at  t  =  0  in  phase  space  will  be  mapped  to  a  set  with 
identical  symplectic  volume  for  any  t  >  0.  More  for¬ 
mally,  Liouville’s  theorem  appears  in  two  equivalent 
formulations,  the  local  and  the  global  form  [23]. 

Theorem  1  (Liouville’s  theorem  (local  form)) 

Let  f*(x)  be  the  flow  of  a  vector  field  f  on  a  chart 
X  of  a  Hamiltonian  system,  i.e.  3 H  :  X  — ■>  M  with 
dim(A’)  =  21V ,  N  £  N  such  that 

f(x)  =(  -  °  liVnXJV")  DxH{x)  Vx£X  (5) 

Then  for  all  x  £  X  and  for  all  times  t  for  which  the 
flow  is  defined, 

Dxffx)  £  Sp2JV;  det  ( Dxf(x ))  =  1  (6) 

(  Sp 2n  denotes  the  group  of  symplectic  matrices  of  size 
2 N  x  2 N.)  The  matrix  of  partial  derivatives  of  the  flow 
with  respect  to  the  initial  conditions  x  is  symplectic 
and  its  determinant  is  one. 

The  global  form  states  that  /*  maps  a  measurable  set 
of  initial  conditions  to  a  set  of  equal  measure. 

Definition  1  (Volume  preservation)  A  map  S  : 

X  — >  X  is  locally  volume  preserving  at  a  point  x  £  X  if 
|  det  ( DxS(x ))  |  =  1.  Its  local  volume  at  x  is  defined  to 
be  det  (DxS(x)).  It  is  volume  preserving  (or  globally 
volume  preserving)  if  |  det  ( DxS(x ))  |  =  1  \/x  £  X . 

This  definition  retains  the  familiar  informal  notion 
of  volume  preservation  (the  “global”  integral  version 
arises,  after  all,  from  the  local  determinant  condition) 
at  the  expense  of  a  slight  degree  of  imprecision  in 
terminology.  Upon  cursory  inspection,  it  might  be 
thought  that  conservative  “piecewise  holonomic”  [25] 
systems  automatically  satisfy  the  hypotheses  of  Liou¬ 
ville’s  theorem.  By  fixing  f  at  a  particular  but  ar¬ 
bitrary  time  t,  a  “degenerate”  hybrid  dynamical  sys¬ 
tem  can  be  defined  on  a  single  chart  X\  =  X  with 
one  vectorfield  f\  =  f  and  the  threshold  function 
hi(ft(xo),xo,t)  =  t  —  t.  The  resulting  stride  map 
S  =  F\  =  f(l(-)  with  t\  =  t  then  obviously  satisfies 
det  (DXS( x))  =  1  Vx  £  X.  However,  for  a  threshold 


equation  that  is  not  purely  time-dependent  but  also 
depends  on  /*( xo )  and  Xq,  the  evolution  time  t\  is  de¬ 
pendent  upon  the  initial  condition:  t±  =  t\( xo),  and 
det  (Dxf^ixo))  ±  1  in  general.  Hence  for  a  gen¬ 
eral  hybrid  dynamical  system  in  which  the  threshold 
functions  are  not  purely  time-dependent,  the  determi¬ 
nant  of  the  Jacobian  of  the  stride  map  S  (4)  cannot 
be  expected  to  be  of  absolute  value  one,  even  if  all  the 
vector  fields  are  Hamiltonian  and  all  transition  func¬ 
tions  are  volume  preserving. 

Liouville’s  theorem  precludes  the  asymptotic  stabil¬ 
ity  of  a  Hamiltonian  system,  since  an  asymptoti¬ 
cally  stable  equilibrium  point  reduces  a  finite  phase 
space  volume  to  a  single  point.  This  would  require 
lim^oo  det  (Dxft(x))  =  0  for  all  x  in  the  basin  of 
attraction  of  the  asymptotically  stable  equilibrium 
point.  However,  because  Liouville’s  theorem  is  not 
guaranteed  to  apply,  asymptotic  stability  of  piecewise- 
defined  holonomically  constrained  conservative  Hamil¬ 
tonian  systems  whose  discrete  time  behavior  can  be 
described  by  an  appropriate  projection  of  a  stride  map 

5.3  has  been  observed  in  the  literature.  Examples  in¬ 
clude  a  discrete  version  of  the  Chaplygin  sleigh  [25,  26] 
and  low-dimensional  models  of  legged  locomotion  in 
the  horizontal  and  sagittal  planes  [11,  13,  14].  In  all 
of  those  cases,  some  threshold  functions  are  not  solely 
time-dependent  and  the  stride  map  is  not  volume  pre¬ 
serving  -  a  necessary  condition  for  asymptotic  stabil¬ 
ity.  In  particular,  at  an  asymptotically  stable  fixed 
point  x,  |  det(Dx5(S))|  <  1. 

Having  established  the  non-applicability  of  Liou¬ 
ville’s  theorem  to  general  hybrid  dynamical  systems, 
we  will  present  criteria  in  the  next  section  under 
which,  nevertheless,  the  volume  preservation  property, 

|  det  (DxS(x))  |  =  1  does  indeed  hold.  The  result  could 
be  called  a  point  Liouville  theorem  for  stride  map  fixed 
points,  because  in  distinction  to  the  local  form  of  Liou¬ 
ville’s  theorem,  which  holds  for  all  points  of  symplectic 
phase  space,  our  theorem  only  holds  at  fixed  points, 
x,  of  S. 

2.3  A  point  Liouville  theorem  for  hy¬ 
brid  dynamical  systems 

In  order  to  prove  that  |  det  (DxS(x))  \  =  1  at  a  fixed 
point  of  S,  additional  assumptions  and  an  additional 
structure  of  the  underlying  vectorfields  fa  is  needed. 
In  particular,  we  require  that  the  vectorfields  fa  pos¬ 
sess  a  time  reversal  symmetry  (for  a  survey  of  time 
reversal  symmetries  in  dynamical  systems  see  [27];  for 
an  extensive  review  see  [28]): 

3The  term  “piecewise  holonomic  system”  was  introduced  in 

[25]. 
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Definition  2  (Time  reversal  symmetry)  A  vec¬ 
tor  field,  f  on  a  chart  X  admits  a  time  reversal  sym¬ 
metry  G  :  X  — >  X  with  G  an  involution1  (G  o  G  =  id) 

if 

DxG-f=-foG.  (7) 

or,  equivalently ,  if 

Go  f*  =  f~*  o  G  .  (8) 

We  will  next  introduce  a  further  property  of  the  stride 
map  factors,  Sa,  of  S'  =  S2  o  Si,  namely  that  they  can 
be  written  as  time  reversed  flow  maps  Sa  =  Ga  o  Fa 
or  Sa  =  Fa  o  Ga.  We  restrict  our  investigation  to  a 
subset  of  fixed  points  of  S,  namely  the  ones  that  are 
also  fixed  points  of  the  time  reversed  flow  maps  Sa, 
Such  fixed  points  we  will  call  symmetric  in  analogy 
to  certain  fixed  points  of  reversible  diffeomorphisms 
(see  definition  6  in  Appendix  C.l).  Fixed  points  of 
this  kind  will  be  shown  to  lie  on  distinguished  orbits 
termed  symmetric  [29];  Such  orbits  have  been  recog¬ 
nized  in  the  prior  legged  locomotion  literature  as  use¬ 
ful  steady  state  target  trajectories  in  the  control  of 
one-legged  hoppers  [9]  and  also  serve  as  steady  state 
target  trajectories  in  this  paper. 

Definition  3  (Symmetric  orbit  of  a  time  re¬ 
versible  vector  field) 

The  orbit  of  a  vector  field  f  with  time  reversal  sym¬ 
metry  G  is  called  symmetric  if  it  is  invariant  under 
G  [29].  This  definition  of  symmetric  orbits  coincides 
with  the  notion  of  neutral  orbits  introduced  in  [9]  and 
formalized  in  [30]. 

Theorem  2  Let  x  be  a  fixed  point  of  Sa  =  Ga  o  Fa, 
where  Fa  is  the  flow  map  of  a  vector  field  fa  with  time 
reversal  symmetry  Ga.  Then  x  lies  on  a  symmetric 
orbit  of  fa. 

Proof:  If  x  is  a  fixed  point  of  Sa  then  there  exists 
a  time  t  such  that  Ga  o  ff(x)  =  x.  If  a:  lies  on  a 
symmetric  orbit  then  Vi  G  [0,f|  3 1'  G  [0,i]  :  f]  (x)  _= 
Ga  °  fa(x).  Let  t'  =  t-t.  Then  f£(x)  =  Ga  o  /*“*  o 
Ga(x)  =  Ga  ofl(x). 

Clearly,  S  locally  preserves  volume  at  a  symmetric 
fixed  point  x  if  its  time  reversed  flow  maps  do.  On 
the  other  hand,  involutions  are  known  to  be  volume 
preserving  at  their  fixed  points: 

Theorem  3  The  determinant  of  the  Jacobian  of  an 
involution  G  :  Rw  — >  at  a  fixed  point  x  of  G  is 

4In  this  paper,  we  restrict  ourselves  to  involutive  time  rever¬ 
sal  symmetries,  although  a  more  general  definition  can  be  found 
in  [27], 


plus  or  minus  one. 

Proof: 

GoG  =  id 

D(G  o  G)  (x)  =  lNxN  \/x£Rn 
DG(G(x))  ■  DG(x)  =  lNxN  (9) 

Since  G(x)  =  x,  (9)  implies  that: 

DG(x)-DG(x)  =  ljvxiv 
=>  clet2(DG(5))  =  1  (10) 

Hence  a  criterion  for  Sa  being  an  involution  is  needed. 

Lemma  1  If  ta  is  Sa  invariant,  that  is,  ta  o  Sa  =  ta 
on  a  set  Xha ,  then  Sa  o  Fa  is  an  involution  on  Xha  ■ 
Proof:  Let  Xq  G  Xha. 

Sa  o  Sa(xo)  — 

Ga  o  Fa  o  Ga  o  Fa(x 0) 

Gao/Ms«W)oGao/^N(xo)  =  (11) 

/-t0(Sa(xo))o/t„(x0)(a.o)  =  a.0. 

A  condition  for  the  Sa  invariance  of  ta  is  now  given, 
in  turn,  as  follows. 

Lemma  2  A  necessary  condition  for  the  Sa  invari¬ 
ance  ofta  is  ha(Ga(x0),Ga  o  ftaa(-Xo\x0),ta(x0))  =  0 
Vico  G  Xha . 

Proof:  If  ta  is  Sa-invariant  then  ta{xf)  must  solve  the 
threshold  equation  for  Sa(x 0): 

haif^  O  Gq  O  /*<»(* »)(®0),Ga  °  ftiXo)(x0),ta(x0))  = 

ha(Ga(x0),Ga  O  /M*°>(®o),ia(so))  =  0  (12) 

Assuming  that  ta{ Xq)  is  also  the  minimal  solution  of 
the  threshold  equation  for  Sa(x 0),  it  follows  that  the 
condition  of  Lemma  2  is  also  sufficient,  and  one  con¬ 
cludes  that  t, a  is  invariant  under  Sa.  Lemma  2  es¬ 
sentially  checks  that  the  threshold  function  ha  “pre¬ 
serves”  the  time  reversal  symmetry  of  fa. 

Combining  Lemmas  1  and  2  and  Theorem  3  consti¬ 
tutes  a  point-wise  Liouville  theorem  for  discrete  sys¬ 
tems  of  the  form  S  =  S20S1  at  symmetric  fixed  points. 
The  generalization  to  a  stride  map  composed  of  more 
than  two  time  reversed  flow  maps  Sa  is  straightfor¬ 
ward.  As  a  final  observation  that  we  will  require  be¬ 
low  (in  Appendix  A),  note  that  if  Theorem  3  has  been 
shown  to  hold  for  Sa  =  Ga  o  Fa,  it  also  holds  for 
reverse  time  flow  maps  of  the  form  Sa  =  Fa  o  Ga : 

Lemma  3  If  Sa  =  Ga  o  Fa  is  an  involution,  then 
S' a  =  Fao  Ga  is  an  involution,  too. 
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Figure  1:  Coordinate  convention  of  SLIP  with  pitch¬ 
ing  dynamics.  In  the  text,  the  COM  coordinates 
will  be  parametrized  by  cartesian  coordinates,  i.  e. 
y  =  Csin(^>)  and  z  =  ((cos (ip).  In  flight,  the  leg  an¬ 
gle  (j>  will  in  general  be  a  function  of  time  and  of  the 
SLIP’s  liftoff  state:  </>(t,x o). 

Proof: 

S' a  o  S' a  =  FaoGaoFao  Ga 

=  (Ga  o  Ga)  oFao  Ga  o  Fao  Ga 

" - v - ' 

—id 

=  Ga  o  Ga  =  id  . 

2.4  SLIP  dynamics 

2.4.1  Modeling  assumptions 


tion  of  the  COM  in  terms  of  cartesian  (y,  z )  and  po¬ 
lar  ((  =  \J y2  +  z2 ,  ip  =  arctan (y/z))  coordinates  with 
the  coordinate  origin  at  the  foothold.  The  body  is  as¬ 
sumed  to  remain  in  the  sagittal  plane  and  its  config¬ 
uration  is  parametrized  by  SE( 2)  coordinates  (y,  z,  9) 
or 

Trajectories:  A  full  stride  consists  of  a  stance  and  a 
flight  phase:  in  stance,  we  assume  the  foothold  is  fixed, 
the  leg  compressed  and  the  body  moves  in  the  positive 
y  direction  y  >  0;  in  flight,  the  body  describes  a  ballis¬ 
tic  trajectory  under  the  sole  influence  of  gravity.  The 
stance  phase  starts  with  the  leg  uncompressed  and 
ends  when  the  leg  has  reached  its  rest  length  (  again. 
Then  the  flight  phase  begins  and  ends  when  the  mass¬ 
less  leg  -  appropriately  placed  -  touches  the  ground. 
Stability  investigations  in  this  paper  are  confined  to 
trajectories  that  are  in  the  vicinity  of  symmetric  tra¬ 
jectories  in  both  stance  and  flight,  where  e.g.  the  liftoff 
and  touchdown  vertical  heights  are  equal. 

Control:  No  continuous  control  is  exerted  during 
stance  and  flight,  the  corresponding  vector  fields  do 
not  change  from  stride  to  stride.  The  only  con¬ 
trol  authority  consists  in  determining  the  transitions 
between  flight  and  stance  by  specifying  the  stance 
and  flight  time.  The  stance  time  is  implicitly  deter¬ 
mined  by  requiring  the  leg  to  undergo  a  compression- 
decompression  cycle,  hence  the  only  designable  control 
authority  consists  in  specifying  the  flight  time,  which 
can  be  implicitly  parametrized  by  the  free  leg  angle 
trajectory  <j)(t,x o)  and  hence  touchdown  COM  height. 
Due  to  the  massless  assumption,  the  leg  can  be  arbi¬ 
trarily  placed  during  flight  at  no  energetic  cost. 


In  this  section  we  establish  the  specifics  of  the  SLIP 
models  considered  in  this  paper.  They  are  listed  in 
terms  of  the  categories:  geometry,  potential  forces, 
control,  and  orbits: 


Potential  forces: 

PI  The  potential  energy  is  given  by  Ep  =  z  + 
V(y,z,6). 


Geometry:  The  3  DOF  sagittal  plane  SLIP  model 
is  shown  in  Fig.  1.  It  shows  a  rigid  body  of  mass  fh 
and  moment  of  inertia  /  with  a  massless  springy  leg 
with  rest  length  (q  attached  at  a  hip  joint  that  coin¬ 
cides  with  the  center  of  mass  (COM).  The  strength  of 
gravity  is  g.  The  approximation  of  a  leg  with  zero  mass 
avoids  impact  losses  at  touchdown  and  simplifies  the 
control.  For  convenience,  all  of  the  following  expres¬ 
sions  are  formulated  in  dimensionless  quantities,  i.e. 

2_  r  —  J-  r  -  — I—,  0  =  6, 


~’Z=Z>Z  = 


t  =  t\  f,y=f,y  = 

V  <.0  <.0  V4 09  ■»*>  V<.0 9 

and  I  = 

g  mC,l 

9  with  respect  to  the  horizontal  and  the  parametriza- 


9  = 


Also  shown  are  the  pitch  angle 


P2  The  non-gravitational  potential  V  is  analytic 
and  satisfies  the  symmetry  relation  V(y,  z,  9)  = 
V(—y,z,—9).  This  condition  does  not  seem  to 
severely  restrict  our  choice  of  potentials,  and  it 
includes  the  often-used  radial  spring  potential 
V(y,z,9)  =  K(C)  f°r  the  2  DOF  model. 

P3  V  factorizes  as  V(y,z,9)  =  Vr(QVp(y,  z,9)  with 
14(1)  =  0.  This  ensures  that  V  is  zero  at  touch¬ 
down  and  liftoff.  Because  of  the  masslessness  of 
the  leg,  V  stays  zero  during  flight. 

After  having  listed  SLIP’s  modeling  assumptions,  we 
will  define  the  stance  and  flight  components  of  the 
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hybrid  SLIP  system  and  identify  time  reversal  sym¬ 
metries  present  in  its  vectorfields. 

2.4.2  Definition  of  the  hybrid  SLIP  system 

The  SLIP  system  consists  of  two  phases  -  stance  and 
flight  -  hence  X  =  {1,2}  with  1  referring  to  stance  and 
2  referring  to  flight.  In  both  phases,  we  choose  the 
same  parametrization  of  the  configuration  space  -  by 
the  cartesian  coordinates  of  the  mass  center  relative 
to  the  fixed  toe,  y,  z,  and  the  orientation  in  of  the 
body  in  the  inertial  frame,  9.  Hence,  both  charts  are 
equal,  Aj  =  X2  =  I?2  x  S1  x  R3  =:  X  with  phase  space 
elements  denoted  by  x  =  (y,  z ,  9 ,  y,  z,  9)T . 

Stance  The  stance  vector  field  reads 

V 
z 

9 

—dyV(y,  z,  9) 

-1  -  dzV (y,  z,  9) 

-)d9V( [y,z,0) 

With  P2  this  vector  field  is  also  analytic  in  x  and 
hence  its  flow  f\(x)  is  analytic  in  t  and  x.  Using  P3 
/i  admits  the  linear  time  reversal  symmetry 

G\  =  diag(— 1, 1,  —1, 1,  —1, 1)  .  (14) 

(the  linear  time  reversal  symmetry  of  (13)  without 
pitching  dynamics  was  already  recognized  in  [30]). 
With  the  “radius”  function  :  x  i— >  \Jy2  +  z2,  the 
threshold  function  is  given  by 

hi(x(t),x0,t)  =  C(x(t))  -  C(x0)  (15) 

Flight  The  flight  vector  field  reads 

h(x)  =  (y,  z,  9, 0,  -1,  o)  (16) 

whose  analytic  flow  is  trivially  computed  as 

/  2/o  +  yot  \ 

z0  +  z0t  -  U 

flix  0)=  9o  +  9ot  (17) 

2/o 

io  - 1 

\  90  J 

Solving  eq.  (7)  with  fa,  the  diagonal  linear  involutive 
time  reversing  symmetry  G2  of  (16)  is  not  uniquely 
defined  and  is  given  by 

G$  =  diag(Tl,  1,  Tl,  ±1,  -1,  ±1)  •  (18) 


As  will  become  clear  later  in  the  next  section,  in  or¬ 
der  to  define  a  stride  map  as  in  (4),  the  time  reversal 
symmetries  should  match  for  stance  and  flight,  hence 
G 2  =  G\  =:  G  is  chosen. 

The  threshold  function  h2  for  a  general  leg  placement 
parametrized  by  the  angular  trajectory  xq)  (see 
Fig.  1)  becomes  zero  when  the  toe  touches  the  ground 

h2(x(t),x0,t)  =  z(t)  -  cos{(/>(t,x0))  ■  (19) 

and  implicitly  defines  the  control  input  t2{x o).  If  4>  de¬ 
pends  on  aio  -  the  liftoff  coordinates,  feedback  control 
is  employed.  The  design  of  the  function  <p  constitutes 
the  control  authority  in  our  SLIP  model. 

2.4.3  Discrete  time  behavior  of  SLIP  locomo¬ 
tion:  Poincare  section,  return  map,  and 
controlled  plant  model 

Poincare  section  A  SLIP  stride  consists  of  stance 
and  flight,  therefore  its  stride  map  should  be  written 
as  S  =  F2o  F\.  The  end  of  the  stance  phase  is  charac¬ 
terized  by  the  liftoff  event,  detected  by  the  threshold 
equation  hi ;  the  end  of  flight  is  characterized  by  the 
touchdown  event,  detected  by  the  threshold  equation 
h2 .  The  factorization  of  S  suggests  a  Poincare  section 
V  that  is  the  surface  of  the  touchdown  event,  where 
the  leg  length  is  one  and  the  COM  is  to  the  left  of  the 
foothold: 

V  =  {x  £  X  :y2  +  z2 -l  =  0,y  <0}  .  (20) 

Return  map  We  would  like  to  factor  S  into  time 
reversed  flow  maps  Sa  in  order  to  satisfy  a  prerequisite 
of  Lemma  1.  This  is  accomplished  by  inserting  the 
square  of  the  common  time  reversal  symmetry  G: 

S  =  F2oGoGoFi  (21) 

However,  S  does  not  formally  constitute  a  return  map 
for  the  Poincare  section  V,  because  as  detailed  in  Sec¬ 
tion  2.4.1,  trajectories  of  relevance  to  forward  locomo¬ 
tion  have  a  monotonically  increasing  fore-aft  compo¬ 
nent;  y(t),  hence,  cannot  be  periodic.  On  the  other 
hand,  there  is  an  effective  projection  informally  built 
into  the  SLIP  modeling  assumption  P3.  At  the  be¬ 
ginning  of  stance,  the  y-coordinate  of  the  coordinate 
origin  must  be  reset  to  the  new  foothold  in  order  to 
interpret  Vr  as  a  radial  leg  potential  (or,  more  awk¬ 
wardly,  one  could  reset  the  definition  of  the  potential 
function  at  each  new  touchdown).  Both  issues  can  be 
resolved  by  projecting  out  the  y-entry  of  S.  A  fur¬ 
ther  dimensional  reduction  is  possible  because  of  con¬ 
servation  of  energy  in  both  stance  and  flight  phase. 


Formally,  the  total  energy 

E(x(t))  =  \(y2(t)  +  z2(t)  +  ie2(t))  + 

z(t)  +  V(y(t),z(t),0(t)) 

=:  E0 

can  be  interpreted  as  a  constant  parameter  of  the  SLIP 
system  and  can  then  be  used  to  eliminate  the  y  vari¬ 
able  y(t)  =  E~{t)  (Eq),5  with  x  being  the  projection  of 

x  onto  its  “non -y,y”  components:  II :  X  — >  X\  x  i— > 
x  =  (z,9,  z,9)T .  A  return  map  R  acting  on  the  re¬ 
duced  Poincare  section  X  =  M  x  S1  x  I2  with  inde¬ 
pendent  coordinates  x  can  then  be  written  as 


Using  a  leg  angular  trajectory  to  implement  feedback 
control,  the  threshold  equation  implicitly  defines  the 
flight  time  ^(/s)  by 

h(k)  =  min {t :  h-2{fl{G  o  R1{x{k))),  y(k),  t)  =  0} 

(26) 

Using  the  explicit  form  of  /12,  (19),  this  expression  for 
the  flight  time,  in  turn,  is  a  function  of  the  control 
input 

u{k)  =  H{y{k))=f>{-,y{k))  (27) 

where  H  parametrizes  the  leg  angle  trajectory  in  terms 
of  the  output  vector  y(k)  and  the  “dummy”  variable 
t,  denoted  by  •. 


The  y  and  y  components  of  F-2  and  G  are  completely 
decoupled  from  the  other  components,  hence  the  pro¬ 
jector  II  can  be  pulled  to  the  right  in  order  to  define 
two  return  map  factors  Ra 


fl  =  F2oGonoGoFl0S  , 


where  F2  and  G  are  the  obvious  restrictions  of  F2  and 
G  to  the  reduced  Poincare  section  X .  If  Sa  are  in¬ 
volutions,  we  want  the  involutive  character  to  per¬ 
sist  for  Ra.  This  is  obvious  for  R2  =  S2 ■  For  R\ 
it  requires  E  o  II  =  id  on  the  range  of  G  o  F\  o  E. 
Let  x\  —Go  F\  o  E(cco)  with  xo  €  V.  y\  is  the  G- 
reflected  ^-coordinate  at  liftoff,  hence  yi  =  —  \/l  —  zf; 
and  yi  =  E~^(E0).  Therefore  j/i  =  Eo  II(yi)  and  R\ 
is  an  involution  if  S\  is  one. 


Controlled  plant  model  Having  defined  the  closed 
loop  return  map  on  the  reduced  Poincare  section,  we 
clarify  the  relation  of  this  closed  loop  return  map  to 
the  controlled  plant  model  formalism  introduced  in 
Section  1.2.  Since  the  control  parameter  of  our  SLIP 
model  is  the  flight  time  and  quantities  used  for  feed¬ 
back  are  the  liftoff  coordinates,  the  controlled  plant 
model,  introduced  conceptually  above  (2),  can  now  be 
written  in  touchdown  coordinates  as 

x(k  + 1)  =  fl2(k)  oGoR^xik)) 

y(k)  =  C(G  o  R^xik)))  (25) 

5Given  an  equation  g(y,x)  =  go,  the  corresponding  implicit 
function  will  be  written  as  y  =  g^ 1  (go )  ■ 


2.4.4  Notation 


The  salient  symbols  used  in  this  paper  are  next  listed, 
with  brief  explanations  of  their  meanings. 


General  hybrid  system  definitions 

1 

finite  index  set,  enumerated  by  a 

Xa 

chart:  phase  space  of  a  dynamical  system 

t,  X 

time,  chart  element  (dimensionless) 

'U 

vector  field  of  a  dynamical  system  on  Xa 

fi 

flow  of  fa  on  Xa 

Fa 

flow  map 

~TrP 

a 

transition  function 

ha 

threshold  function:  triggers  chart  transition 

ta(x  o) 

evolution  time  on  chart  Xa  starting  at  Jo 

V 

Poincare  section  (surface  in  Xa ) 

Xa 

reduced  Poincare  section 

Ra 

return  map  factor  on  Xa 

R 

return,  Poincare  map 

In  general,  an  element  or  a  map  without  the  diacritic 
T  denotes  an  element  of  the  reduced  Poincare  section 
Xa  or  a  map  on  Xa. 


Other  definitions 

Ga 

involutive  time  reversal  symmetry 

xha 

set  where  partial  stride  map  is  an  involution 

sa 

stride  map  factor  on  Xa 

s 

stride  map 

n 

projector  from  Xa  to  Xa 

E 

map  from  Xa  to  Xa 

V 

conservative  SLIP  potential  without  gravity 

3  Stability  and  control  of  SLIP 
models 

In  this  section  the  stability  and  control  of  SLIP  models 
will  be  analyzed  via  the  return  map  R  and  its  factors 
Ra.  In  Section  3.1  it  is  first  shown  that  the  stance 
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factor  Ri  is  locally  volume  preserving  at  a  fixed  point 
x,  independent  of  the  specific  form  of  the  potential  V 
as  long  as  the  conditions  PI  -  P3  are  satisfied.  We  will 
then  derive  an  expression  for  the  local  volume  of  R2 
as  a  function  of  the  leg  angle  trajectory  </>.  Combining 
these  two  results  will  give  a  necessary  condition  for 
stability  of  a  SLIP  model  in  terms  of  the  controlled  leg 
angle  trajectory  <f>.  Note  that  by  different  SLIP  models 
we  mean  SLIP  models  that  have  potentials  satisfying 
the  conditions  PI  -  P3  but  that  differ  in  their  leg 
angle  trajectories  4>. 

In  the  remaining  portions  of  this  section,  we  use  the 
preceding  analysis  to  explore  an  informal  relation  be¬ 
tween  the  “degree  of  stability”  as  manifest  in  the  sin¬ 
gularity  of  the  linearized  discrete  return  map  and  the 
“cost  of  feedback.”  The  latter  is  judged  with  respect 
to  a  number  of  quantitative  and  qualitative  features 
of  known  relevance  in  robotic  implementations.  These 
informal  “cost”  measures  are  introduced  and  moti¬ 
vated  in  Section  3.2  and  are  shown  to  be  quantifi¬ 
able  using  the  preceding  analysis.  Next,  in  Section 
3.3  we  apply  the  results  of  3.2  to  the  study  of  several 
2  DOF  SLIP  models  (i.e.,  SLIP  models  without  pitch¬ 
ing  dynamics)  that  have  appeared  in  the  literature, 
classifying  them  with  respect  to  the  “cost”  properties 
previously  introduced.  Finally,  in  Section  3.4  we  in¬ 
troduce  a  new  3  DOF  SLIP  model  that  offers  a  more 
realistic  description  of  the  physical  robot  RHex  oper¬ 
ating  under  the  influence  of  its  open  loop  gait  gener¬ 
ating  “clock”  [1],  We  apply  the  analytical  methods  of 
Section  3.1,  characterizing  sensory  “cost”  and  control 
benefit  laid  out  in  Section  3.2,  and  are  able  to  give  for 
the  first  time  conditions  on  the  RHex  clock  param¬ 
eters  -  some  necessary  for  gait  stability,  and  others 
sufficient  for  gait  instability. 

3.1  Computation  of  the  local  return 
map  volume 

3.1.1  Stance 

In  this  section  we  will  apply  the  results  of  Section  2.3 
to  show  that  R\  is  an  involution  by  showing  that  S\  is 
an  involution  for  a  SLIP  model  satisfying  the  assump¬ 
tions  of  Section  2.4.1.  We  first  apply  Lemma  2:  Given 
£ 1  =  ii (a?o),  the  threshold  equation  in  Lemma  2  reads 

hi(G(x0),G(x(ti)),ti)  =  C(£o)  -  C(z(h))  =  0  (28) 

However,  since  this  is  just  the  negative  of  the  origi¬ 
nal  threshold  equation  hi(x(ti),Xo,ti)  =  C(^(H))  — 
£(a?o)  =  0,  H  is  a  solution  of  (28).  Assuming  that 
£1(2:0)  is  indeed  the  minimal  solution  of  the  threshold 
equation  for  5i(xq)  for  all  xq  £  A,  Lemma  2  can  be 


applied  to  prove  that  Si  is  an  involution  on  A^  =  A. 
By  the  arguments  in  Section  2.4.3,  R\  is  also  an  invo¬ 
lution  and  Theorem  3  now  implies  that  R\  is  locally 
volume  preserving  at  its  fixed  point. 

3.1.2  Flight 

We  now  derive  a  formula  for  the  determinant  of  the 
Jacobian  of  the  flow  map  F2  given  an  arbitrary  leg 
angle  trajectory  <j>(t,x 0).  This  is  used  to  compute  the 
determinant  of  the  Jacobian  of  the  partial  return  map 
R2  =  F2  0  G  at  a  fixed  point  of  i?2- 
Note,  in  contrast  to  i?i,  that  |  det(DxR2(x))\  can  be 
computed  directly  for  any  specific  leg  angular  trajec¬ 
tory  <j>  using  the  closed  form  expression  of  the  flight 
phase  flow  (17).  Nevertheless,  in  Appendix  A  Lemma 
2  is  applied  to  a  particular  family  of  leg  angle  trajec¬ 
tories  in  order  to  classify  which  of  the  resulting  flight 
phase  return  maps  are  involutions. 

The  threshold  function  h2  for  a  general  leg  angle  tra¬ 
jectory  (f>  is  h2(x(t),xo,t)  =  z(t)  —  cos(c/>(t, xo))  (19). 
Setting  h2  =  0  determines  the  time  from  leg  liftoff 
(tLO  =  0)  to  leg  touchdown  txD  =  t2.  Because  h2 
is  a  transcendental  map,  a  closed  form  expression  for 
£2(2:0)  cannot  be  found  in  general. 

It  should  be  pointed  out  that  the  dependence  of 
</>(£,  Xq)  on  the  flight  time  £  is  redundant  in  the  sense 
that  the  leg  angle  is  irrelevant  to  the  dynamics  of  the 
system  except  at  the  touchdown  time  £td(2:o)-  Specif¬ 
ically,  a  given  flight  time  £td(£o)  =  £2(2:0)  can  be 
enforced  by  a  purely  state  dependent  leg  angle  “tra¬ 
jectory”  (j>(x 0)  =  arccos  (^(£2(2:0)))  or  by  any  time  de¬ 
pendent  trajectory  <j>'(t,x 0)  that  satisfies 

<^(£2(2:0), x0)  =  0(xo)  .  (29) 

The  advantage  of  including  time  as  an  additional  ar¬ 
gument  of  <f>  will  be  pointed  out  in  Section  3.3.1. 

The  flow  map  F2  takes  the  state  vector  Xo  from  its 
value  at  leg  liftoff  to  that  at  touchdown:  F2(a;o)  = 
x(trD )•  A  fixed  point  of  a  symmetric  flight  trajectory 
satisfies  x  =  S2(x)  =  F2  o  G(x). 

The  determinant  of  the  Jacobian  of  F2(xq)  = 
f2TD<'x°\xo)  can  easily  be  computed  from  the  expres¬ 
sion  for  the  flight  phase  flow  (17),  bearing  in  mind 
that  the  flight  time  tTD( xq)  also  depends  on  the  ini¬ 
tial  conditions: 

det(Dx(F2)(2;o))  = 

1  —  dz0tTD( Xo)  +  ZodZQtTD(Xo)  +  ^O^o^-D^o)  (30) 

This  expression  exemplifies  the  remarks  in  Section  2.2, 
since  it  will  reduce  to  one,  in  general,  only  if  Itd  is 
independent  of  the  initial  conditions  Xq-  Hence  using 
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implicit  differentiation  of  (19)  the  determinant  can  be 
written  in  terms  of  partial  derivatives  of  cf>(t,x o): 

A  1  num 

det  {DxF2(x o))  =  1  +  2lden  (31) 

2  t—tj'D 

with 

A2num  =  sm(<j)(t,x0))  ■ 

{dio(j){t,x 0)  -  z0dZo(l)(t,xo)  -  90d9o(l)(t,xo )) 

+  t  —  Zo 

A2 den  =  sin  x0))dt<j)(t,  x0)-t  +  z0  . 

Albeit  tx d  cannot  be  computed  in  closed  form  in  gen¬ 
eral  because  of  the  transcendental  nature  of  h2,  we 
know  that  at  a  fixed  point  x  of  F2  o  G  with  i’o  := 
G(x)  the  liftoff  and  touchdown  heights  are  identical 
and  hence  txD  =  2zo-  Therefore,  sm(cf)(tTD,  So))  = 

—  \/l  —  Zq  and  6{txD )  =  —  &o-  The  eigenvalues  of  the 
partial  return  map  F2  o  G  at  such  a  fixed  point  are 
{1, 1,  — 1,  —  det(Dx(F2  o  G(S)))}. 

Because  G  =  diag(l,  — 1,  — 1, 1),  the  determinant  of 
the  Jacobian  of  R2  and  F2  are  related  as 

det{DxR2{x))  =  det  {DxF2(G(x)))  (32) 

2  DOF  SLIP  model  For  the  2  DOF  SLIP  model 
without  pitching  dynamics,  the  0,  9  variables  are  ab¬ 
sent  and  F2,  G,  and  R2  are  2-dimensional  maps.  The 
determinant  of  the  flight  phase  flow  map  simplifies  to 

clet(DxF2(x0))  =  1  +  (33) 

sin  (</>(f,  x0))  (dio(t>(t,  xp)  -  zpdZo</>(t,  x0))  +  t-  z0 

sin (</>(t,  x0))dt<f)(t,  x0)-t  +  z0  t=tTD 

The  eigenvalues  of  the  partial  return  map  F2  o  G  at 
its  fixed  point  x  are  {1,  —  det(Dx(F2  oG(i)))}.  With 
G  =  diag(l,  —1),  the  determinants  of  the  Jacobians  of 
R2  and  F2  are  related  as 

det(DxiZ2(a;))  =  -  det  (DXF2  (£(*)))  (34) 

3.1.3  Local  volume  of  the  return  map  at  a 
symmetric  fixed  point 

Having  derived  expressions  for  |  det(Dxi?i(x))|  and 
|  det(DxJ?2(a:))|  in  the  two  previous  sections  at  fixed 
points  x  of  R\  and  R2,  the  composition  of  R  of  those 
two  partial  return  maps  R  =  R2  o  R1  can  be  used  to 
factor  the  determinant  |  det(DR(x))  |  at  a  symmetric 
fixed  point  x,  i.e.  a  fixed  point  that  is  common  to 
both  R\  and  R2  (see  Section  2.3): 

|  det(DxJ?(x))|  =  \det(DxR2(R1(x)))\\det(DxR1(x))\ 

S - - - v - / 

=  1 

=  |det(FxF2(i))|  (35) 


Hence  a  necessary  condition  for  local  asymptotic  sta¬ 
bility  of  R  at  x  is  |  det(Di?(x))|  <  1,  whereas  a 
sufficient  condition  for  local  asymptotic  instability  is 
|  det(-Di?(a;))|  >  l.6  The  factor  |  det(Dxf?2  (x))|  is 
governed  by  the  time  of  flight  (30)  which  in  turn  de¬ 
pends  upon  the  functional  form  of  the  leg  angle  trajec¬ 
tory  <f>  (31).  Demanding  stability  of  R  at  a  symmetric 
fixed  point  therefore  imposes  conditions  on  </>,  or,  us¬ 
ing  the  formalism  of  controlled  plant  models,  on  H oC 
specified  in  (27). 

3.2  Deadbeat  control  and  singular  re¬ 
turn  map  Jacobians 

3.2.1  Control  and  sensor  modeling 

For  discrete  systems,  three  different  degrees  of  local 
stability  can  be  distinguished,  which  are  character¬ 
ized  by  the  eigenvalues  of  the  Jacobian  of  the  closed 
loop  return  map  at  a  fixed  point:  i)  all  eigenvalues 
are  within  the  unit  circle;  ii)  all  eigenvalues  are  within 
the  unit  circle  and  some  are  zero  (“singular  control”); 
iii)  all  eigenvalues  are  zero  (“deadbeat  control”).  In 
general,  the  more  singular  the  closed  loop  return  map 
the  quicker  the  transient  behavior '  but  the  higher  the 
“cost”  of  control  and  the  more  vulnerable  to  model¬ 
ing  errors.  Although  we  are  not  interested  in  pursuing 
formal  optimality  conditions,  assessing  the  overall  sen¬ 
sory  cost  of  various  control  alternatives  is  of  central 
concern  in  physical  robotics  applications.  One  rea¬ 
sonable  approach  that  we  adopt  here  is  to  count  the 
number  and  characterize  the  “quality”  of  the  sensed 
variables  required  to  complete  the  feedback  loop  of 
the  controlled  plant  model  (2).  Here,  “quality”  refers 
to  the  frame  of  reference  of  the  feedback  variables, 
since  body  frame  sensing  is  generally  easier  to  accom¬ 
plish  than  inertial  sensing.  In  SLIP  models,  feedback 
control  is  parametrized  by  the  leg  angle  trajectory 
(f>(t,  Xo),  where  Xo  are  the  state  variables  taken  at  a 
certain  event.  Intuitively,  three  different  aspects  of 
sensory  cost  can  be  readily  distinguished: 

SI  Detection  of  the  event  where  the  feedback  vari¬ 
ables  are  taken 

i)  easy  for  liftoff:  can  be  implemented  in  a  SLIP 
hopper  by  a  simple  switch  at  the  toe 

6 Note  that  necessary  and  sufficient  conditions  for  stability 
would  require  the  knowledge  of  the  eigenvalues  of  R  at  x.  How¬ 
ever,  eigenvalues  of  a  composition  of  two  maps  do  not  factorize 
into  eigenvalues  of  the  two  individual  maps  unless  the  maps 
commute  -  i.e.,  both  are  diagonalizable  via  the  same  similarity 
transformation. 

'  This  is  motivated  by  the  fact  that  a  function  from  P ;V  to 
P :V  whose  Jacobian  has  rank  K  <  N  everywhere  maps  an  N- 
dimensional  volume  to  a  /f-dimensional  volume. 
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ii)  difficult  for  flight  phase  apex:  requires  mea¬ 
surement  of  vertical  velocity  z 

52  Enforcement  of  the  angle  trajectory  <f>:  a  leg  angle 
trajectory  <f>  specified  with  respect  to  an  inertial 
frame  requires  inertial  sensing  for  enforcement 
(i.e.  feedback  control),  as  opposed  to  a  leg  an¬ 
gle  trajectory  specified  with  respect  to  the  body 
frame.8 

53  Sensing  of  the  feedback  variable  xq  by  the  output 
map  C  (2): 

i)  dimension  of  the  domain  (number  of  argu¬ 
ments)  of  C 

ii)  position  versus  velocity  measurement:  posi¬ 
tions  are  in  general  easier  to  measure  than  ve¬ 
locities 

iii)  “quality”:  inertial  versus  non-inertial  (body 
frame)  quantities 

Because  we  exploit  in  this  paper  the  factorization  of 
R  into  stance  and  flight  phase,  it  is  natural  to  work  in 
“liftoff  coordinates”  -  i.e.,  on  the  Poincare  section  V  - 
hence,  the  feedback  variables  are  naturally  assumed  to 
be  taken  at  the  “easily  detected”  liftoff  event  as  noted 
in  SI.  We  appraise  in  Section  3.3.1  the  alternative 
choice  of  working  formally  in  apex  coordinates  (not 
to  be  confused  with  the  physically  unattractive  choice 
of  taking  the  sensory  feedback  measurements  at  the 
apex  event).  Criteria  S2  and  S3  can  be  addressed  by 
rewriting  the  leg  angular  trajectory  <j>  that  is  defined 
in  an  inertial  frame  (see  Fig.  1)  as 

x0)  =  C(x0))  -  9{t)  .  (36) 


where  <f>Bo  is  the  leg  liftoff  angle  with  respect  to  the 
body  normal  (see  Fig.  1)  and  if) bo  is  the  leg’s  angular 
velocity  at  liftoff  measured  in  the  body  frame.  Speci¬ 
fying  this  trajectory  in  the  body  frame  yields 

<j>(t,x0)  =  ^Cb^^BoAbo)  ~  6(t)  (38) 

In  summary,  the  3  DOF  SLIP  model  allows  the  dis¬ 
tinction  of  the  “quality”  of  sensing  required  for  a  par¬ 
ticular  control  input  which  in  turn  enables  an  assess¬ 
ment  of  the  “cost”  of  control. 

3.2.2  Deadbeat  control  requires  singular  re¬ 
turn  map  Jacobians 

In  this  section,  we  observe  that  deadbeat  control  of  a 
2  or  3  DOF  SLIP  model  requires  the  return  map  Ja¬ 
cobian  to  be  globally  singular  -  not  just  at  the  control 
target  fixed  point  x  but  on  the  entire  reduced  Poincare 
section  X. 

For  the  full  nonlinear  closed  loop  plant  model  the  re¬ 
turn  map  R  is  deadbeat  if  there  exists  aA'eN  such 
that 

Rk{x)  =  x\/xeX  (39) 

for  a  specified  target  x.  This  means  that 

DxR(Rk~\x))  •  DxR(Rk~2(x))  .... 

* DxR(x )  =  xdim(A’)  ^  %  (40) 

A  necessary  condition  for  this  is 
det  ( DxR(Rk~'(x )))  =  0  for  some  1  <  i  <  K. 
Since  eq.  (40)  must  be  valid  \/x  €  X ,  we  need 
det  ( DxR(x ))  =0  \/x  €  X. 


The  second  term  in  (36)  indicates  that  <f>c  is  specified 
with  respect  to  the  SLIP’s  body  frame,  as  will  be  the 
case  in  all  3  DOF  SLIP  models  in  this  paper.  For  2 
DOF  SLIP  models,  9  is  not  defined  and  this  term  is 
absent. 

It  is  not  possible  to  distinguish  S3iii,  “quality”  (i.e. 
inertial  vs  non-inertial  frame  based)  in  the  2  DOF 
setting,  since  by  its  very  geometry,  body  frame  co¬ 
ordinates  cannot  be  introduced.  On  the  other  hand, 
the  additional  body  pitch  degree  of  freedom  of  the  3 
DOF  SLIP  model  allows  this  distinction  to  be  made. 
A  leg  angle  trajectory  that  only  uses  sensing  with  re¬ 
spect  to  the  body  reference  frame  S3,  can  be  modeled 
by  the  following  output  map  Cb  ■ 


(arccos(,2o)  +  Oq'' 
V x~zl  / 


Cb{x  o) 


(37) 


®Note  that  this  feedback  control  cannot  be  modeled  straight- 
forwardly  in  our  simplified  SLIP  system  because  of  the  mass¬ 
lessness  of  the  leg. 


3.2.3  General  solution  of  leg  angle  trajectory 
with  singular  return  map  Jacobians 

As  will  be  reviewed  in  Section  3.3.1,  2  DOF  SLIP  mod¬ 
els  with  globally  singular  return  map  Jacobians  have 
featured  prominently  in  the  literature  -  both  dead¬ 
beat  and  non-deadbeat.  In  this  section  we  will  derive 
the  general  form  of  leg  angle  trajectories  that  render 
the  return  map  Jacobian  globally  singular.  In  general, 
the  matrix  DXF2  will  have  full  rank.  If,  under  the  in¬ 
fluence  of  a  particular  leg  angle  trajectory,  </>(f,; To), 
the  second  factor  of  the  closed  loop  return  map  is 
rank  deficient  for  all  state  vectors,  det  (.Da;  #2(20))  = 
0  =  det(DxF2(xo)),  and  if  a  stable  fixed  point  exists, 
then,  as  discussed  in  Section  3.2.1,  one  would  expect 
a  “more  rapid”  convergence  to  this  fixed  point  than  if 
the  matrix  had  full  rank.  Since  formula  (31)  is  valid 
for  arbitrary  flight  times,  not  just  at  a  fixed  point  of 
i?2,  a  partial  differential  equation  for  globally  singu¬ 
lar  leg  angle  trajectories  <j)(t,x 0)  can  be  obtained  by 
setting  (31)  to  zero. 
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First,  setting  eq.  (30)  to  zero  yields  a  partial  differ¬ 
ential  equation  for  the  touchdown  time  txD,  which  by 
the  method  of  characteristics  [31]  has  the  general  so¬ 
lution 

t-TD{zo,  zo,  90, 0o)  =tA  +  t(za,  0a,  9 a)  (41) 

where  r  is  an  arbitrary  differentiable  function  of  the 
three  arguments  and  the  transformation  from  liftoff 
coordinates  to  apex  coordinates  is  given  by 


tA 

=  Z0 

eA 

=  90 

•2 

ZA 

Zr\ 

=  z0+- 

eA 

=  6q  +  9qZo 

(42) 

specifying  the  time  from  liftoff  to  apex,  the  pitch¬ 
ing  velocity  at  apex,  the  apex  height,  and  the  apex 
pitch  angle.  Since  a  scalar  function  tTD  is  deter¬ 
mined  by  this  partial  differential  equation,  the  rank 
of  Dx  F2  will  generally  be  reduced  by  only  one  using 
tA  +  t{za,  9a,  9a)  as  the  flight  time. 

Now  recall  that  the  time  of  flight  function,  trn  =  ti 
arises  in  our  application  as  the  solution  of  an  im¬ 
plicit  function  defined  by  the  leg  angle  trajectory  p 
(19).  Thus,  we  must  next  impose  a  corresponding 
“singularity”  condition  on  (f>  that  guarantees  the  de¬ 
sired  property  in  txD-  Setting  eq.  (31)  to  zero, 
det(DF2(a:o))  =  0,  yields  the  partial  differential  equa¬ 
tion  for  globally  singular  leg  angle  trajectories, 

dt4>{t,xo)+dz0(t>(t,xo)-zodZo(t>{t,xo)-dod80(t)(t,xo)  =  0 

The  general  solution  of  this  linear,  homogeneous,  first 
order  partial  differential  equation  by  the  methods  of 
characteristics  is 

4>(t,x0)  =  ®(t-tA,zA,9A,QA)  (43) 

where  4>  is  an  arbitrary  differentiable  function  of  its 
four  arguments. 

3.3  2  DOF  SLIP  models:  sensor  re¬ 

quirements  and  stability 

This  section  focuses  on  2  DOF  SLIP  models  with  re¬ 
spect  to  sensor  requirements  in  their  feedback  loop. 
First,  it  is  shown  that  all  2  DOF  SLIP  models  with 
globally  singular  return  map  Jacobians  require  a  mea¬ 
surement  of  the  vertical  velocity,  either  explicitly 
through  the  arguments  of  (f>  or  implicitly.  Then  the 
dimensional  reduction  of  the  return  map  that  follows 
from  the  globally  singular  return  map  Jacobians  is  il¬ 
lustrated  with  four  different  2  DOF  SLIP  models  that 


have  already  appeared  in  the  literature.  A  stable  2 
DOF  SLIP  model  with  full  rank  return  map  Jaco¬ 
bian  is  also  presented  to  illustrate  the  power  of  our 
analysis  in  the  low  dimensional  setting.  Since  the  re¬ 
duced  Poincare  section,  X,  is  only  two  dimensional 
for  the  2  DOF  model,  the  presence  of  complex  con¬ 
jugate  eigenvalues  of  the  linearized  return  map  at  a 
given  fixed  point  strengthens  our  stability  criteria  to 
the  point  that  the  determinant  magnitude  condition  is 
both  necessary  and  sufficient  for  asymptotic  stability. 
Thus,  as  we  demonstrate,  by  varying  one  parameter, 
asymptotically  stable,  neutrally  stable,  and  unstable 
behavior  can  be  exactly  assigned. 

3.3.1  All  singular  2  DOF  SLIP  models  require 
velocity  sensing 

In  this  section  several  previously  proposed  [13,  14,  9, 
32,  33]  2  DOF  SLIP  control  strategies  are  reviewed 
with  emphasis  on  their  globally  singular  return  map 
Jacobians.  The  general  solution  for  a  globally  singu¬ 
lar  leg  angle  trajectory  for  the  2  DOF  SLIP  model  is 
obtained  from  (43)  by  omitting  the  pitch  coordinates, 
hence  <f>(t,x o)  =  &(t  —  tA,zA).  But  both  control  input 
arguments  require  the  vertical  velocity  measurement 
zo  when  expressed  in  liftoff  coordinates  (42),  which 
leaves  the  constant  trajectory  Xq)  =  const  as  the 
only  globally  singular  leg  angle  trajectory  without  ex¬ 
plicit  velocity  sensing.  We  will  review  four  2  DOF 
SLIP  models  with  globally  singular  leg  angle  trajec¬ 
tories,  pointing  out  that  even  the  leg  angle  trajectory 
(f>(t,  Xq)  —  const  requires  velocity  sensing  for  its  im¬ 
plementation  as  highlighted  in  criterion  S2. 

Constant  leg  touchdown  angle  policy  The  con¬ 
stant  leg  touchdown  angle  policy  proposed  in  [13,  14, 
34,  32]  has  the  simple  form 

4>(t,  xo)  =  2tt  —  P  :  t  >  tA  (44) 

where  /?  is  a  constant  angle  for  all  strides.  No  sensing 
of  the  feedback  variables  S3  is  required,  hence  the  out¬ 
put  map  C  can  be  taken  to  be  a  constant.  Since  the  re¬ 
turn  map  Jacobian  of  this  SLIP  model  is  globally  sin¬ 
gular,  the  return  map  is  effectively  one-dimensional. 
In  [13]  this  one-dimensional  variable  was  taken  to  be 
the  apex  height,  whereas  in  [14]  the  angle  of  the  touch¬ 
down  velocity  was  chosen. 

9 A  similar  leg  angular  trajectory  for  a  3  DOF  SLIP  model 
was  shown  in  [14]  to  yield  asymptotically  stable  behavior  for 
certain  parameter  values.  Although  not  presented  here,  the 
return  map  factorization  introduced  in  this  paper  can  be  applied 
to  this  model  also  to  show  that  its  stance  phase  is  locally  volume 
preserving  at  a  symmetric  fixed  point  whereas  its  flight  phase 
has  a  globally  singular  return  map. 
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A  Poincare  section  volume  and  the  embedded  one¬ 
dimensional  return  map  domain  is  plotted  in  Fig. 
2a),  where  the  return  map  image  X1  :=  R(X)  with 
X  =  [0.8, 0.99]  x  [—1.5,  —0.1]  is  depicted  by  solid  points 
joined  by  a  black  line.  The  color  of  the  points  matches 
the  color  of  the  inverse  images  f?_1  ( XI(zAi ))  of  these 
points.  The  color  corresponds  to  a  parametrization  of 
the  return  map  image  in  terms  of  the  resulting  apex 
heights  ZAt  ■  Since  a  constant  leg  touchdown  angle 
is  prescribed,  the  touchdown  height  is  constant  and 
the  return  map  image  is  a  vertical  line  in  (zo,zo)- 
coordinates.  The  curved  black  line  denotes  the  one¬ 
dimensional  manifold  of  all  possible  fixed  points  for 
arbitrary  leg  angle  trajectories.10  Although  (f>  is  a  con¬ 
stant  and  does  not  explicitly  depend  on  the  velocity 
measurement  of  Zo,  vertical  velocity  sensing  is  implicit 
in  the  derivation  of  the  return  maps  in  [13,  14,  34,  32], 
because  the  leg  angle  is  not  held  constant  through¬ 
out  the  flight  phase,  but  is  assumed  to  be  set  to 
2tt—/3  in  a  time  interval  ( io —  \J  z^  +  2(zo  —  sin  (3),  io  + 
\J Zq  T  2(zo  —  sin/3))  in  which  the  COM  is  above  the 
touchdown  height  sin  / 3 .  Before  this  time  interval  is 
reached,  the  leg  is  assumed  to  be  at  an  angle  where  it 
does  not  interfere  with  the  ground. 

Raibert  controller  The  leg  placement  strategy 
proposed  by  Raibert  [9]  for  a  two  degree  of  freedom 
SLIP  reads 

x0)  =  2t r  -  arcsin  +  ky(y0  -  y)j  (45) 

where  ts  is  the  duration  of  the  stance  phase,  ky  is 
a  feedback  gain  and  y  is  the  desired  forward  speed. 
In  Raibert’s  physical  implementations,  the  duration 
of  the  current  stance  phase  was  approximated  by  the 
measured  duration  of  the  previous  stance  phase.  Here, 
we  will  consider  ts  a  constant.  In  (45)  the  average 
forward  stance  speed  used  in  [9]  was  approximated  by 
?/o-  Now  yo  can  be  expressed  as  yo  =  y/2(E  —  za). 
Hence  eq.  (45)  is  of  the  form  (43)  and  the  return 
map  domain  is  a  one-dimensional  manifold  which  is 

10 By  Theorem  2  a  fixed  point  of  the  time  reversed  stance 
flow  map  5*1  lies  on  a  symmetric  orbit  of  its  vector  field  f\. 
Symmetric  orbits  must  contain  a  fixed  point  of  G  [30]  and  can 
therefore  be  characterized  for  the  2  DOF  SLIP  model  by  the 
two-dimensional  fixed  point  set  FixG  =  {r  G  :  y  =  0,  i  =  0}. 
Fixing  the  energy  Eq  removes  one  dimension,  hence  the  set  of 
all  possible  fixed  points  of  the  return  map  factor  R\  forms  a  one¬ 
dimensional  manifold  in  X.  Given  that  any  x  =  (z,  z)T  with 
i  >  0  lies  on  a  symmetric  orbit  of  the  flight  phase  vector  field 
on  the  reduced  Poincare  section  X,  the  set  of  all  possible  fixed 
points  of  the  return  map  R  is  identical  to  the  one-dimensional 
manifold  of  possible  fixed  points  of  R±.  The  fixed  points  of  R 
are  then  given  by  the  intersection  of  this  line  with  the  return 
map  image. 


depicted  in  Fig.  2b).  The  output  map  for  this  leg 
angular  trajectory  reads  C{x o)  =  za- 

Leg  retraction  and  “optimized  selfstabiliza¬ 
tion”  In  the  leg  retraction  schemes  proposed  in 
[15,  33],  the  leg  is  set  at  a  fixed  angle  a  a  at  the  apex 
of  the  flight  phase  and  then  starts  rotating  towards 
the  ground.  Before  reaching  the  apex,  the  leg  angle 
can  be  arbitrarily  placed  as  long  as  its  toe  does  not 
touch  the  ground.  In  [15],  a  constant  angular  velocity 
to  is  used  (leg  retraction),  i.  e. 

4>{t,x0)  =  aA  +  uj(t-tA)  :t>tA  (46) 

whereas  in  [33]  a  nonlinear  angular  trajectory  that  is 
constant  over  all  strides 

x0)  =  a(t  -tA)  :  t>tA  (47) 

is  employed.  In  both  cases,  the  output  map  is  C(x o)  = 
t-A-  Clearly,  these  two  leg  placement  schemes  are  also 
of  the  form  (43)  and  therefore  the  return  map  im¬ 
age  is  a  one-dimensional  manifold.  These  return  map 
images  are  plotted  in  Figs.  2c  and  cl  respectively. 
Both  return  maps  converge  to  the  same  point,  how¬ 
ever,  the  second  trajectory  [33]  achieves  convergence 
to  a  desired  apex  height  within  one  stride.11  Since  the 
apex  Poincare  section  in  [33]  is  only  one-dimensional 
and  one  control  parameter  -  the  touchdown  time  or 
rather  the  leg  touchdown  angle  -  is  available,  the  de¬ 
sired  apex  height  can  be  reached  within  one  stride. 
On  the  other  hand,  the  touchdown  Poincare  section 
parametrized  by  (z,  z)  is  two-dimensional  and  dead¬ 
beat  control  can  only  be  achieved  within  at  least  two 
strides.  This  seems  to  be  a  contradiction,  since  the 
discrete  time  behavior  of  identical  physical  systems 
parametrized  by  different  Poincare  sections  must  be 
conjugate,  i.e.  related  by  a  coordinate  transforma¬ 
tion.  Particularly  the  dimension  of  the  return  maps 
of  both  parametrizations  must  agree.  In  Appendix  B 
it  is  shown  that  if  all  coordinates  of  the  dynamical 
flow  are  taken  into  account,  the  apex  and  touchdown 
return  maps  are  indeed  conjugate.  However,  because 
the  open  loop  system  is  dynamically  decoupled  in  apex 
coordinates  (i.e.  the  second  variable  does  not  influ¬ 
ence  the  evolution  of  the  first  in  these  coordinates), 
restricting  the  feedback  to  depend  upon  the  first  vari¬ 
able  yields  effectively  a  one-dimensional  closed  loop 
return  map.  This  one-dimensional  nature  is  illus¬ 
trated  in  Fig.  2d),  where  the  one-dimensional  man¬ 
ifold  X1  :=  R(X)  is  plotted  together  with  color-coded 
inverse  images  R_1  ( XI{zAi ))•  As  can  be  seen  in  Fig. 

11  The  angular  trajectory  a  was  obtained  by  numerical  inver¬ 
sion  of  the  apex  height-to-apex  height  return  map  in  order  to 
implement  deadbeat  control. 
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2d,  X1  is  aligned  with  one  of  the  inverse  images,  hence 
in  the  first  stride  an  arbitrary  point  (z,z)  is  mapped 
onto  X1 ,  whereas  in  the  second  stride  all  points  on 
this  manifold  are  mapped  to  the  target  point. 

The  authors  of  [33]  call  this  control  scheme  “optimized 
selfstabilization,”  indicating  a  computational  or  sen¬ 
sory  advantage  over  regular  deadbeat  control.  In  reg¬ 
ular  deadbeat  control,  the  leg  angle  </>  would  be  a  func¬ 
tion  of  both  zq  and  Zo,  requiring  the  sensing  of  both 
liftoff  variables  and  the  online  computation  or  storage 
of  a  lookup  table  for  a  function  from  a  two-dimensional 
to  a  one-dimensional  space.  In  (47)  only  the  sensing 
of  tA  =  Zo  and  a  clock  is  required,  and  a  is  a  function 
from  a  one-dimensional  to  a  one-dimensional  space.  In 
this  context,  “selfstability”  seems  to  refer  to  the  fact 
that  the  leg  angle  is  a  function  of  time  (starting  at 
apex)  only  and  and  does  not  explicitly  depend  upon 
the  liftoff  variables  (zo,  io)i  it  does  not  mean  that  no 
sensing  (e.  g.  detection  of  the  apex)  is  required.  In  the 
next  paragraph  we  address  the  explicit  parametriza- 
tion  of  this  one-dimensional  return  map  manifold  and 
show  how  it  can  be  used  to  reduce  the  sensory  require¬ 
ments  of  control. 


Sensory  requirements  of  globally  singular  con¬ 
trol:  Given  a  globally  singular  2-DOF  SLIP  return 
map  with  leg  angle  trajectory  <j){t,x o),  this  leg  angle 
trajectory  can  be  rewritten  as  Xo)  =  $(t  —  tA,  za) 
according  to  the  results  in  section  3.2.3.  The  corre¬ 
sponding  output  map  can  be  chosen  to  be  C(a;o)  = 
(t_4, za)T ■  This  does  not  constitute  a  sensory  advan¬ 
tage  over  xo  because  still  one  position  and  one  veloc¬ 
ity  measurement  are  required.  The  threshold  function 
reads 

h2(x(t),x0,t)  =  z(t)  -  cos  ($(t  -  tA,zA))  (48) 

=  za-^—a^- - cos  (®(t-tA,zA)) 

Setting  h2  to  zero  implicitly  defines  a  function  A  tA 
with  the  substitution  t  —  tA  — >  A tA(zA).  A tA{zA)  en¬ 
codes  the  direct  control  parameter  during  flight  -  the 
total  flight  time  tA  +  A tA(zA).  A  different  angular 
trajectory  enforcing  the  same  total  flight  time  for  all 
initial  conditions  zq,  io  can  then  be  defined  by  the  in¬ 
verse  A tjj1:  <I>(t  —  tA )  :=  4>(f  —  tA,  A t~^(t  —  t a))  with 
a  new  output  map  C{ xq)  =  tA  whose  only  output  is 
the  flight  time  measured  from  apex.  Hence  a  leg  angle 
trajectory  xo)  that  initially  required  the  sensing 
of  (tA,Z4)T  and  time  can  be  replaced  by  one  that 
only  requires  sensing  of  the  apex,  i.  e.  tA  =  zq,  and 
time.  This  rewriting  of  the  angular  trajectory  makes 
use  of  the  invariance  of  the  flight  time  with  respect  to 
certain  parametrizations  of  <f>  (29)  and  demonstrates 


why  deadbeat  control  for  SLIP  models  can  be  achieved 
with  reduced  feedback  sensing  S3i. 


3.3.2  A  nonsingular,  stable  2  DOF  SLIP 
model  without  velocity  sensing 

We  will  now  investigate  a  2  DOF  SLIP  model  with  a 
full  rank  return  map  Jacobian  where  we  address  both 
S3i  and  S3ii  in  that  no  velocity  sensing  is  required 
for  the  feedback  loop.  For  certain  parameter  values, 
this  model  does  exhibit  asymptotic  stability.  In  the 
previous  2  DOF  examples  of  Section  3.3.1,  once  sin¬ 
gularity  has  been  imposed,  the  determinant  of  the  re¬ 
turn  map  Jacobian  vanishes  and  the  factor  analysis 
can  contribute  no  more  information  to  the  stabiliza¬ 
tion  problem.  However,  as  this  example  shows,  since 
the  return  map  has  dimension  two,  if  we  operate  in  a 
regime  where  the  eigenvalues  are  known  to  have  non¬ 
zero  imaginary  components,  then  the  properties  of  the 
determinant  completely  determine  stability.  We  can 
then  dictate  the  stability  properties  through  a  closed 
form  expression  and  this  is  indeed  how  the  present 
example  has  been  adjusted. 

The  leg  angle  trajectory  for  this  model  reads 

c i>{t ,  Xq)  =  ut  +  k  arccos(20)  +  ola  (49) 


where  u>,  k,  and  a  a  are  constants.  Note  that  z0  does 
not  appear  in  (49),  hence  the  output  map  could  be 
written  as  C(xq)  =  Zq.  For  k  =  1  and  a  a  =  0,  the  leg 
rotates  clockwise  at  a  constant  rate  ui  starting  with  the 
liftoff  angle  arccos(zo)-  This  can  be  considered  a  crude 

2  DOF  SLIP  version  of  the  leg  angle  profile  specified 
by  RHex’s  open  loop  controller  [1],  A  more  elaborate 

3  DOF  SLIP  version  of  RHex’s  open  loop  controller 
will  be  presented  in  Section  3.4.2.  Using  (33)  the  de¬ 
terminant  of  the  Jacobian  of  R  at  a  symmetric  fixed 
point  becomes: 


|  det(Dxi?(x))| 


|clet(D,F2(G(x)))| 
—z(k  —  1) 


(50) 
<  1 


In  order  to  illustrate  the  predictive  power  of 

(50),  we  numerically  approximate  the  determinant 

det(DxR(x))  of  the  full  return  map  for  fixed  SLIP 

parameters  En  =  E°~  =  2.1,  7  =  13,  and  fixed  re- 

mg(  0 

circulation  parameters  ola  =  7r,  u>  =  14  for  different 
k  G  {1/6, 0.5, 1,  2, 3.3}.  Here,  E0  is  the  dimensionless 
total  conserved  energy  of  the  system  and  the  dimen¬ 
sionless  spring  potential  is  U(C)  =  (7/ 2)(£  —  l)2.  We 
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Figure  2:  One-dimensional  return  map  domains  and  their  inverse  images  for  rank-deficient  SLIP-controllers:  a) 
Fixed  leg  angle  touchdown,  b)  Raibert,  c)  Leg  retraction,  d)  Two-step  deadbeat.  All  elements  of  a  colored  line 
in  the  (^o,  io)-plane  are  mapped  to  the  point  with  identical  color.  The  union  of  all  these  points  constitutes  the 
return  map  image.  The  color  corresponds  to  a  parametrization  of  the  return  map  image  in  terms  of  the  resulting 
apex  heights  Za,  ■  The  range  of  apex  heights  considered  is  Za  €  [0.92,1.8].  The  curved  black  line  identical  in  all 
four  figures  denotes  the  set  of  all  possible  fixed  points,  as  explained  in  Footnote  10. 


then  compare  these  values  to  the  values  of  the  deter¬ 
minant  obtained  by  inserting  the  numerically  deter¬ 
mined  fixed  points  x  =  (z  z)T  into  (50).  The  de¬ 
terminants  obtained  in  those  two  different  ways  are 
plotted  in  Fig.  3a  and  agree  to  a  high  precision 
(||  det(Da;i?(x))|  -  |det(DxF2(G(x)))||  <  10"7).  In 
Figs.  3b-d  iterations  of  the  return  map  in  (zq  io)-space 
are  shown  for  k  £  {1/6, 1,3.3}  and  initial  conditions 
off  the  fixed  point.  The  eigenvalues  are  complex  con¬ 
jugate  pairs  in  all  three  cases,  hence,  the  magnitude  of 
the  eigenvalues  computed  in  (50)  specifies  sufficient  as 
well  as  necessary  conditions  for  stability  and  instabil¬ 
ity  in  this  case.  For  k  =  3.3,  formula  (50)  specifies  an 
unstable  fixed  point,  and,  indeed,  the  plot  of  a  numer¬ 
ical  simulation  in  Fig.  3b  depicts  a  typical  trajectory 
spiraling  away  from  a  small  neighborhood  as  required. 


For  k  =  1/6,  formula  (50)  specifies  asymptotic  stabil¬ 
ity,  and  trajectories  spiral  towards  the  fixed  point,  as 
depicted  in  Fig.  3c.  For  k  =  1,  formula  (50)  suggests 
neutral  stability  and  numerical  simulation  verifies  that 
all  trajectories  lie  on  deformed  circles  around  the  fixed 
point  as  plotted  in  Fig.  3d. 

Fig.  3d  is  reminiscent  of  KAM-tori  of  area-preserving 
2D  mappings  (see  [35]),  however,  as  can  be  seen  in  Fig. 
4,  the  phase  space  volume  is  not  preserved  away  from 
the  fixed  point  x  for  k  —  1.  In  Appendix  C  we  invoke 
reversibility  [36]  in  place  of  area-preservation  to  show 
that  the  numerically  observed  neutral  stability  for  the 
leg  recirculation  scheme  with  k  =  1  is  expected. 
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Figure  3:  a)  Comparison  of  the  numerically  computed  determinant  |  det(Da.R(x))|  (+)  of  the  return  map  Jacobian  to  to  the 
determinant  |  det(Da:i?2(G(a;)))|  (o)  obtained  by  using  the  numerically  determined  fixed  points  in  (50).  b)-d)  Trajectories  around  a 
fixed  point.  Because  of  slow  convergence,  only  every  9th  iteration  in  plot  b)  and  every  5th  iteration  in  plot  c)  is  shown. 


3.4  3  DOF  SLIP  models:  body-frame 

sensing  and  stability 


In  this  section  the  control  possibilities  with  body 
frame  sensing  are  explored  for  3  DOF  SLIP  models. 
First,  the  unique  3  DOF  SLIP  model  with  the  body 
frame  sensor  model  (37)  and  a  globally  singular  return 
map  Jacobian  is  presented.  By  comparing  the  number 
of  available  design  parameters  of  this  SLIP  model  to 
the  dimension  of  the  reduced  Poincare  section,  dead¬ 
beat  control  is  excluded.  Then  a  nonsingular  3  DOF 
SLIP  model  with  only  body  frame  sensing  is  intro¬ 
duced.  It  is  modeled  within  the  limitations  of  the  3 
DOF  SLIP  dynamics  after  the  open-loop  controller  of 
RHex  [1]  and  is  shown  to  have  asymptotically  stable 
operating  regions.  A  necessary  condition  for  the  sta¬ 
bility  of  this  model  in  terms  of  a  RHex  clock  parameter 
is  derived. 


3.4.1  Body  frame  sensing  does  not  admit 
deadbeat  control 

We  want  to  investigate  the  possibility  of  deadbeat 
control  with  a  leg  angle  trajectory  of  the  form  (38) 
o)  =  </>cB(L<iW^Bo)-0(i)>  i-e-  using  only  body 
fame  sensor  information  in  the  feedback  loop  and  spec¬ 
ifying  the  leg  angle  trajectory  in  the  body  frame. 

As  shown  in  Section  3.2.2,  deadbeat  control  requires 
globally  singular  return  map  Jacobians  and  hence 
—  0(t)  must  be  cast  into  the  form 
$(t-tA,ZA,0A,dA)  (43).  While  9(t)  =  dA  +  dA(t-tA) 
does  satisfy  this  functional  form,  (f>cB  (L  ^b0>  ^b0) 
does  not,  except  for  (j>cB  ( t,4>B0,4>B0 )  =  const.  We  will 
present  numerical  evidence  in  the  form  of  an  asymp¬ 
totically  stable  trajectory  at  particular  parameter  val¬ 
ues  of  a  3  DOF  SLIP  model  in  order  to  show  that 
stable  behavior  is  possible  with  the  leg  angle  trajec¬ 
tory 

<p(t,  xo )  =  27 t  —  /3  —  9(t)  .  (51) 

This  3  DOF  SLIP  model  bears  close  resemblance  to 
the  LLS  model  [11]  of  horizontal  legged  locomotion, 
where  at  the  end  of  each  stance  phase  the  new  stance 
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Figure  4:  Contour  plot  of  |  det(Dxii(:r))|  for  the  leg 
angle  trajectory  (49)  with  k  =  1.  The  fixed  point 
x  «  (0.8772,  — 0.0764)t  lies  on  the  |  det(.DxiZ(®))|  =  1 
contour,  whereas  volume  (area)  is  not  preserved  away 
from  the  fixed  point. 


leg  is  set  at  a  fixed  angle  with  respect  to  the  non- 
inertial  body  axis,  thus  implementing  a  similar  leg  an¬ 
gular  trajectory.12  A  sample  discrete  trajectory  on  the 
3-dimensional  Poincare  section  is  shown  in  Fig.  5  for 
a  potential  of  the  form  P3  with  Vj-(£)  =  (y/2)(^  —  l)2 
and  Vp(ip,  9)  =  1  +  cgg92  +  cg^O^  +  c^ip2.  The  mo¬ 
tivation  for  this  potential  is  discussed  in  Section  4.1. 
Given  that  the  only  design  parameter  of  (51)  is  (3,  a 
target  point  ( z ,  0,  9)  in  the  reduced  three-dimensional 
Poincare  space  cannot  be  specified  a  priori.  Hence  the 
possibility  of  deadbeat  control  for  the  3  DOF  SLIP 
model  with  the  body  frame  sensor  model  (37)  must 
be  discarded. 


3.4.2  A  RHex-like  3  DOF  SLIP  model  with 
body  frame  sensing 

In  this  section  a  leg  placement  strategy  for  the  con¬ 
trol  of  the  3  DOF  SLIP  model  with  full  rank  is  in¬ 
vestigated.  Its  importance  lies  in  the  fact  that  this 
leg  placement  strategy  is  modeled  after  the  open-loop 
controller  employed  in  RHex  [1]  within  the  limitations 
of  the  3  DOF  SLIP  model.  The  angular  reference  tra¬ 
jectories  prescribed  by  RHex’s  open  loop  clock  con¬ 
troller  [1]  are  specified  by  the  (dimensionless)  param¬ 
eters  tc,ts,<pa,<po',  see  Fig-  6-  For  one  half  of  the  clock 
period  tc,  the  trajectories  for  the  left  (L)  and  the  right 
(R)  tripod  can  be  expressed  as  functions  of  time  in  the 


12Note,  however,  that  the  flight  duration  is  zero. 


Figure  5:  A  sample  discrete  trajectory  on  the 
3-dimensional  Poincare  section  parametrized  by 
(zo,6o,9o)  converging  to  an  asymptotically  stable 
fixed  point.  Because  of  the  rank-deficient  nature  of  the 
leg  placement,  z o  is  a  function  of  the  other  Poincare 
section  variables  Zo  =  cos(9q  +  7t/2  —  (3).  The  values 
of  the  dimensionless  variables  characterizing  this  sys¬ 
tem  are  egg  =  400,  eg ^  =  —12,  =  0,  E$  =  2.1, 

7  =  13.25,  I  =  0.489,  and  (3  =  1.0562. 


robot’s  body  frame  (B)  as 


■At)  = 


LOst  +  ip0 


0  <  t  < 


Uft+W-%-)  +  <Po 


(52) 


l(t)  = 


•  +  UJ  ft  +  Ifio 


Uf 


0  <t< 

tc-t 


2  — 


<t  <  f 


where  us  =  <  uif  =  .  These  angular  trajec¬ 

tories  are  depicted  in  Fig.  6.  They  are  enforced  at 
each  leg  of  the  robot  by  a  simple  PD-controller. 

In  order  to  implement  this  controller  in  a  simplified 
3  degree  of  freedom  SLIP  controller,  the  following  as¬ 
sumptions  are  made: 


R1  RHex’s  clock  should  prescribe  motions  with  (sub¬ 
stantial)  flight  phases,  i.e.  ts  <  tc/2. 

R2  During  stance  the  (virtual)  stance  leg  can  be  ap¬ 
proximated  by  a  SLIP;  this  means  that  there  is 
no  net  torque  apart  from  gravity  on  the  leg. 

R3  The  PD-controller  that  enforces  the  angular  ref¬ 
erence  trajectories  (52)  during  flight  has  infinite 
gains  and  tracks  those  trajectories  without  errors. 
In  consequence,  as  the  present  stance  leg  lifts  off, 
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4>br  (t)  =  4>br  (t)  -  <po  4>bl  ( t )  =  4>bl  ( t )  -  <fo 


with  respect  to  the  inertial  frame  is  given  by 


Figure  6:  Illustration  of  liftoff  and  touchdown  events 
in  the  body  frame  of  a  RHex-inspired  leg  recirculation 
scheme  [1]  for  SLIP  with  pitching  dynamics. 

the  PD-controller  can  be  assumed  to  have  posi¬ 
tioned  the  second  (present  flight)  leg  at  the  exact 
angle  with  respect  to  the  first  (present  stance)  leg 
as  specified  by  (52).  In  contrast,  during  stance, 
the  angular  position  and  velocity  variables  evolve 
according  to  the  SLIP  stance  mechanics  (13). 

Assumption  R1  simply  focuses  attention  on  RHex’s 
dynamical  regime  as  opposed  to  the  possibly  quasi¬ 
static  operation  available  to  platforms  with  sufficiently 
high  leg  number.  Assumption  R2  is  justified  by  ex¬ 
periments  and  simulation  studies  of  RHex  operating 
in  the  relevant  dynamical  regime  [37,  2,  34].  Assump¬ 
tions  R2  and  R3  make  the  controlled  SLIP  model  a 
pseudo-clock  controlled  system,  where  the  clock  signal 
is  turned  off  during  stance  and  turned  on  at  a  reset 
phase  at  liftoff.  In  the  following  derivation,  the  left 
leg  (L)  is  chosen  to  be  the  stance  leg. 

In  order  to  cast  the  reference  trajectories  in  an  ex¬ 
pression  for  (/>(t,  Xq)  using  assumption  R3  we  need  to 
express  (j>BR  (the  flight  leg  angle  with  respect  to  the 
body  frame)  in  terms  of  the  angle  of  the  stance  leg  at 
liftoff  4>BL{tLO )-13  In  addition,  we  need  to  transform 
leg  angles  specified  in  the  body  frame  to  the  respective 
angles  in  the  inertial  frame.  The  relation  between  the 
body  and  the  inertial  frame  is  given  by 

4>b  =  4>  +  0  (53) 

Hence  the  leg  angle  trajectory  of  the  right  (flight)  leg 

13Note  that  the  time  in  c/>bl  differs  from  the  time  in  <f>(t,  lo); 
it  is  not  reset  at  the  beginning  of  a  new  flight  phase. 


<t>R{t)  =  4>BR{t  +  tLO)  -  6{t)  =  4>(t,  x0)  (54) 

where  <lo  =  +  $o)  is  the  time  with  respect 

to  the  RHex  clock  when  liftoff  occurs  and  4>l0  = 
arccos(zo)  is  the  angle  of  the  stance  leg  at  liftoff  in  the 
inertial  frame.  This  procedure  is  illustrated  within  the 
body  frame  in  Figure  6. 

Using  the  expressions  for  the  RHex  clock  trajectory 
(52),  we  obtain  an  angular  trajectory  for  the  3  DOF 
SLIP  system  that  at  any  instant  in  time  has  the  gen¬ 
eral  form 

(j)(t,  cco)  =  u>t  +  k  (arccos(2o)  +  &o)  +QA  —  0(t)  (55) 

V - V - / 

=0BO 

with  u>  £  {ujf,cus},  and  k  £  { ,  1,  the  dif¬ 
ferent  values  of  can  easily  be  derived  and  are 

not  important  in  this  context.  This  expression  is 
of  the  form  (38)  with  the  body  frame  sensor  model 
Cb[x o)  =  arccos(zo)  +  9q  =  4>b0i  he.  no  body  frame 
velocity  measurement  is  required.  The  functional  form 
of  eq.  (55)  translates  into  different  functional  expres¬ 
sions  for  |  det(DF2(G(ir)))|  which  are  distinguished  by 
the  location  of  liftoff  and  touchdown  with  respect  to 
the  piecewise-linear  leg  angle  trajectories.  We  enu¬ 
merate  these  six  cases  by  two  numbers  (LO  — >  TD) 
which  denote  the  region  in  Fig.  6  where  liftoff  and 
touchdown  occurs. 


(LO  ->  TD) 

|det(£>F2(G(®)))| 

(1-1) 

(1-2) 

.  w-*-2 

(-Z  +  6y/l  —  Z2)-\-CJfy/l  —  Z2 

(2-2) 

(1-3) 

1 

(2-3) 

(3-3) 

1  H - : — ■  - <  1 

(  —  Z+Qy/l  —  Z2  )  +CJS  y/l  —  Z2 

Table  1:  Functional  expressions  for  j  det(T>F2(G(x)))| 
for  different  locations  of  the  liftoff  and  touchdown 
event. 

Based  on  the  properties  of  symmetric  orbits  that  we 
focus  on  in  this  paper,  a  necessary  condition  on  a 
RHex  clock  parameter  for  asymptotic  stability  can 
now  be  derived. 

A  necessary  condition  for  asymptotic  stability 

For  symmetric  orbits,  the  liftoff  and  touchdown  leg 
angles  in  the  inertial  frame  are  of  equal  magnitude 
but  opposite  sign:  ^l0  =  — 4>R(trD )■  This  also  holds 
for  the  pitching  angles:  6l0  =  — 9(tTD )•  Using  (53)  to 
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translate  the  liftoff  and  touchdown  angles  to  the  body 
frame 

<t>BL{tLo)  =  4>lo+0lo 
4>BR(tTD)  =  (Pr^td)  +  d(tTo) 

we  obtain  <pBL{tLo)  +  4>BR(tTD)  =  0.  With  the  def¬ 
initions  (pBL  =  (pBL  ~  V o  and  4>br  =  (f>BR  ~  as  in 
Figure  6,  this  equality  can  be  rewritten  as 

< t>BL(tLO )  +  4>BR{tTD)  =  —2(^0  (56) 

By  locating  the  liftoff  and  touchdown  angles  of  (56) 
for  symmetric  orbits  within  regions  1-3  in  Figure  6,  a 
table  of  possibly  allowed  liftoff-touchdown  transitions 
as  a  function  of  the  sign  of  RHex  clock’s  leg  offset 
angle  tpo  can  be  derived: 


(LO  ->  TD) 

<Po  <  0 

O 

II 

O 

9- 

yo  >  0 

(1-1) 

(1-2) 

no 

no 

yes 

(2-2) 

(1-3) 

yes 

yes 

yes 

(2-3) 

(3-3) 

yes 

no 

no 

If  9  >  ,  z  ,  the  determinant  for  the  transitions 

y/1  —  Zz  7 

(2  — >  3)  and  (3  — ■>  3)  will  be  less  than  one.  Hence  a 
necessary  condition  for  asymptotic  stability  is  ipo  <  0 
provided  that  9  >  On  the  other  hand,  a  suf¬ 

ficient  condition  for  instability  (including  neutral  sta¬ 
bility)  is  (po  >  0  provided  that  9  > 

It  should  be  emphasized  that  the  expressions  in  table  1 
are  independent  of  the  specific  3  DOF  SLIP  potential 
V  as  long  as  the  general  conditions  PI  -  P3  listed  in 
2.4.1  are  obeyed.  However,  a  specific  SLIP  model  does 
influence  the  location  of  fixed  points  x  and  therefore 
the  numerical  value  of  the  determinant  for  the  cases 
1  — >  1, 1  — >  2, 2  — »  3, 3  — »  3  as  well  as  the  location  of 
the  eigenvalues. 

This  analysis  provides  for  the  first  time  a  partial  ex¬ 
planation  for  the  surprising  stability  of  the  open-loop 
clock  driven  robot  RHex. 

4  Application:  Toward  hierar¬ 
chical  control  of  a  hexapedal 
robot 

In  this  section  we  explore  numerically  the  applicabil¬ 
ity  of  these  ideas  to  the  robot  RHex  [1].  In  Section  4.1 
we  review  the  general  approach  to  hierarchical  control 


[5]  in  the  specific  context  of  a  RHex-like  anchor  sys¬ 
tem  compared  to  the  SLIP-template  with  a  physically 
motivated  stance  phase  potential  and  leg  angle  tra¬ 
jectory.  In  Section  4.2  we  present  simulation  results 
and  assess  the  degree  of  correspondence  between  the 
simulated  anchor  and  its  putative  template. 


4.1  Control  of  the  anchor  by  control¬ 
ling  the  template 

A  template  is  a  low  dimensional  model  of  a  mecha¬ 
nism  operating  within  a  specified  environment  that  is 
capable  of  expressing  a  specific  task.  To  anchor  this 
low  dimensional  model  in  a  more  physically  realistic 
higher  degree  of  freedom  representation  of  the  robot 
and  its  environment,  we  seek  controllers  whose  closed 
loops  result  in  a  “prescriptive”  correspondence  (de¬ 
fined  in  Section  1.2)  of  the  dynamics  of  the  high  and 
low  degree  of  freedom  models.  Hence  the  controller 
must  a)  force  the  high-dimensional  anchor  to  follow 
the  dynamics  of  the  low-dimensional  template  (“an¬ 
choring”),  and  b)  control  the  template  to  achieve  a 
certain  task.  In  the  case  at  hand  the  anchor  is  given 
by  the  robot  RHex,  whereas  the  template  is  given  by 
a  3  DOF  SLIP  model  as  in  Section  2.4.1  with  a  leg 
angle  trajectory  (55).  Note  that  this  template-anchor 
hierarchy  includes  the  intrinsic  abstraction  of  neglect¬ 
ing  lateral  dynamics  and  focusing  on  the  sagittal  plane 
motion.  The  analysis  of  the  previous  chapters  will  be 
shown  to  give  insights  into  part  b);  we  assume  that 
the  SLIP-anchoring  mechanism  (part  a))  has  already 
been  addressed  by  either  deliberate  design  [6]  or  by 
the  interaction  of  the  controlled  robot  with  its  envi¬ 
ronment,  as  has  been  shown  for  RHex’s  steady  state 
behavior  [37,  2]. 

As  a  first  step  towards  the  goal  of  devising  a  new  con¬ 
troller  for  RHex  or  optimizing  RHex’s  open  loop  con¬ 
troller,  we  provide  numerical  evidence  for  the  agree¬ 
ment  of  the  stability  properties  of  a  RHex-like  model 
programmed  in  the  SimSect  simulation  environment 
[7]  with  the  stability  properties  of  the  3  DOF  SLIP 
model  introduced  in  Section  3.4.2.  In  particular, 
we  show  that  in  a  physically  interesting  operational 
regime,  stable  simulations  in  SimSect  correspond  to 
stable  fixed  points  of  the  corresponding  SLIP  model. 
Here,  “correspondence”  is  established  by  fitting  sim¬ 
ulation  data  to  the  3  DOF  SLIP  model  with  stance 
phase  potential  introduced  in  Section  3.4.1  and  with 
the  RHex-like  leg  angle  trajectory  (55)  using  the  RHex 
clock  parameters  of  the  SimSect  simulation.  This 
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stance  phase  potential 


4.2.1  Correspondence  of  trajectories 


V(y,z,0)  =  (1  +  ceeO2  ±  cg^dip  +  c^ip2) 

(57) 

is  the  generalization  of  a  potential  V {ip,  9)  =  (7/2)  (£— 
1)2(1  +  c-e+^{9  ±  ip)2)  which  could  be  implemented  in 
a  physical  single-leg  hopping  robot  by  a  passive  leg 
spring  and  a  passive  torsional  spring  between  the  body 
and  the  leg.  The  generalized  potential  (57)  would  re¬ 
quire  passive  springs  attached  to  a  “co-moving”  in¬ 
ertial  frame,  preventing  a  physical  implementation 
with  a  one-legged  robot.  However,  we  believe  that 
this  model  reasonably  approximates  the  moments  due 
to  the  “outrigger”  front  and  back  legs,  compressed 
against  the  horizontal  ground  surface,  in  the  tripod 
stance  phase  of  RHex. 


4.2  Correspondence  between  RHex- 
like  simulations  and  their  fitted 
SLIP  models 

Simulations  of  a  RHex-like  hexapedal  robot  were  run 
in  SimSect  [7]  over  a  discretized  range  of  clock  pa¬ 
rameters  of  RHex’s  open  loop  clock  controller  [1] 
that  respect  assumptions  Rl,  R2,  R3  of  Section 
3.4.2:  tc  G  [0.235, 0.245],  <ps  G  [0.84, 1.04],  ip0  G 
[—0.16, 0.04],  df  G  [0.52,0.6],  where  the  duty  factor  df 
is  defined  as  df  =  (tc  —  ts)/tc.  Of  the  resulting  1815 
SimSect  simulations,  522  (=  29%)  were  stable  accord¬ 
ing  to  the  criteria  of  Appendix  D.l.  Then  a  3  DOF 
SLIP  model  with  the  stance  potential  (57)  and  the 
leg  angle  trajectory  (55)  was  fit  to  those  stable  cases 
following  the  fitting  procedure  outlined  in  Appendix 
D.2. 

As  is  detailed  below,  the  3  DOF  SLIP  model  approx¬ 
imates  the  24  DOF  SimSect  steady  state  dynamics 
surprisingly  well  given  the  gulf  in  dimension.  Specif¬ 
ically,  the  trajectory  fitting  errors  are  very  small  on 
average,  the  fixed  points  of  the  SimSect  simulations 
and  the  fitted  SLIP  models  are  within  the  same  order 
of  magnitude,  and  the  asymptotic  behavior  agrees  in 
almost  all  cases. 

However,  while  the  fitted  SLIP  models  provide  good 
correspondence  once  a  specific  SimSect  operating 
point  has  been  selected,  it  is  not  the  case  that  a  priori 
specification  of  clock  parameters  yields  a  SLIP  model 
whose  fixed  point  locus  and  stability  predicts  that  ob¬ 
served  in  the  SimSect  model.  In  this  sense,  the  present 
SLIP  model  provides  a  descriptive  but  not  prescriptive 
representation  of  the  SimSect  dynamics. 


The  quality  of  the  fit  is  assessed  for  each  stable  sim¬ 
ulation  by  the  two  fitting  error  numbers  A yzL2  and 
A 9l2  as  described  in  Appendix  D.3.  The  average  fit¬ 
ting  error  and  standard  deviation  for  both  errors  for 
the  522  stable  SimSect  simulations  for  the  cartesian 
coordinates  is  small  A yzL2  =  3.82  ±0.42%  and  of  sim¬ 
ilar  magnitude  as  the  fitting  errors  for  2  DOF  SLIP 
models  observed  in  [34],  whereas  the  average  fitting 
error  for  the  pitch  coordinates  A 9l2  =  93.65  ±25.76% 
is  considerably  larger.  As  an  illustration  of  the  fitting 
results  a  sample  SLIP  fit  is  presented  in  Figure  7.  The 
data  trajectories  of  y{t),y{t),z{t),z{t),9{t),9{t),({t) 
are  plotted  together  with  the  trajectories  of  their  fit¬ 
ted  SLIP  models. 
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Figure  7:  Trajectories  of  a  stable  SimSect  simulation 
and  the  corresponding  trajectories  of  the  fitted  SLIP 
model.  The  fitting  errors  are  A yzL2  =  3.56%  and 
A  6L2  =  51.46%. 


The  large  A 9l2  fitting  error  is  an  indication  that  the 
proposed  3  DOF  SLIP  model  is  not  a  sufficiently  accu¬ 
rate  abstraction  of  SimSect’s  pitching  dynamics.  An¬ 
other  contributing  factor  to  the  size  of  A 6l2  is  the 
fact  that  the  magnitudes  of  both  the  9{t)  and  9{t) 
trajectories  are  small,  which  makes  the  denominator 
of  the  fitting  error  A 9l2  (64)  small.  Another  devia¬ 
tion  of  the  dynamics  of  the  fitted  SLIP  model  from 
SimSect  is  apparent  in  the  y{t)  trajectories  of  Fig.  7 
which  is  typical  of  all  fitted  SLIP  models.  Here,  the 
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fitted  SLIP  trajectory  is  out  of  phase  with  the  SimSect 
COM  trajectory;  nevertheless  the  fitting  error  A i/zl2 
is  small  because  of  the  large  average  value  of  the  Sim¬ 
Sect  trajectory  that  enters  (64).  We  believe  that  the 
acceleration  in  the  forward  direction  during  the  leg 
compression  phase  seen  in  SimSect  as  opposed  to  a 
deceleration  in  the  corresponding  SLIP  model  is  due 
to  the  non-conservative  nature  of  the  SimSect  model, 
where  energy  is  pumped  into  the  system  by  the  hip 
torques,  and  then  lost  to  damping  and  friction. 

4.2.2  Correspondence  of  fixed  points 

Next  we  investigate  whether  the  fixed  point  of  the  fit¬ 
ted  SLIP  model  accurately  predicts  the  fixed  point  of 
the  corresponding  SimSect  run.  In  general,  the  fitted 
SLIP  model  with  the  initial  condition  obtained  from 
its  SimSect  simulation  as  described  in  Appendix  D.3 
will  not  operate  at  a  fixed  point  of  its  return  map  R. 
Hence  a  root-finding  algorithm  (Matlab’s  ‘fsolve’)  is 
employed  to  determine  the  fixed  point  Sslip  (if  it  ex¬ 
ists)  of  the  return  map  R  of  the  fitted  SLIP  model  as 
well  as  the  eigenvalues  {Ai}i=x,...,4  of  the  Jacobian  of 
the  return  map:  DxR(xslip)-  The  fixed  point  Sslip  is 
then  compared  to  the  “fixed  point”  Ssim  -  the  appro¬ 
priate  projection  of  the  initial  data  point  of  the  Sim¬ 
Sect  simulation  stance  data.  Scatter  plots  of  the  com¬ 
ponents  of  the  fixed  points  Sslip  and  Ssim  of  all  stable 
SimSect  simulations  are  shown  in  Fig.  8.  For  perfect 
correspondence  of  the  SimSect  and  SLIP  dynamics,  all 
fixed  point  components  should  lie  on  the  identity  lines. 
While  the  orders  of  magnitude  of  the  components  of 
the  fixed  points  match,  the  components  are  in  gen¬ 
eral  not  well  correlated,  except  for  zslip  and  5sim  if  a 
constant  offset  is  taken  into  account.  The  fixed  points 
■ZSim  assume  an  almost  constant  value  and  the  pitching 
components  of  xslip  are  very  close  to  zero  and  under¬ 
estimate  the  magnitude  of  the  pitching  components  of 

^Sim* 


•2SLIP  ^SLIP 

0.87  0.89  -0.2  -0.1 


Figure  8:  Scatter  plots  of  components  of  the  fixed 
points  of  stable  SimSect  simulations  versus  the  cor¬ 
responding  fixed  point  components  of  the  fitted  SLIP 
models.  For  perfect  correspondence,  all  points  should 
lie  on  the  identity  line. 


considered.  On  the  other  hand,  if  the  magnitudes  of 
the  individual  eigenvalues  are  only  slightly  larger  than 
one,  e.g.  |A;|  <  1.05,*  =  1,  ...,4  as  is  the  case  for  92% 
of  the  522  stable  SimSect  runs,  then  the  instability 
might  only  be  revealed  after  many  iterations  (see 
Fig.  3b  for  a  “weakly”  unstable  trajectory  of  a  2 
DOF  SLIP  model  with  |Ai|  =  |A2|  «  1.007).  Hence 
the  predicted  instability  of  the  fitted  SLIP  models 
for  SimSect  simulations  with  certain  RHex  clock 
parameters  might  not  be  discernible  from  stability 
given  the  criteria  in  Appendix  D.l  due  to  the  finite 
amount  of  simulation  time  for  each  simulation  and 
and  also  due  to  the  limitations  of  the  SLIP-SimSect 
correspondence  as  discussed  above. 


4.2.3  Correspondence  of  stability  at  a  fixed 
point 

Given  the  numerically  determined  eigenvalues 
{ Ai}i=i;...j4  of  a  fitted  SLIP  model  at  its  fixed  point 
5slip,  its  local  asymptotic  stability  properties  can 
be  assessed.  The  magnitude  of  the  determinant 
|n|=1Ai|  agrees  to  a  high  numerical  precision  with 
the  appropriate  expression  in  Table  1,  predicting 
instability  for  ip0  >  0  and  allowing  asymptotic 

stability  for  ipo  >  0.  However,  SimSect  simulations 
that  are  stable  in  the  sense  of  Appendix  D.l  are  found 
with  similar  frequency  of  occurrence  for  ipo  >  0  as 
well  as  for  ipo  <  0  for  the  range  of  clock  parameters 


5  Conclusions 

In  this  paper  we  use  the  example  of  the  SLIP  locomo¬ 
tion  model  to  show  how  factored  analysis  of  the  return 
map  may  be  a  useful  new  tool  in  the  stability  analysis 
of  hybrid  Hamiltonian  systems.  Specifically,  we  derive 
a  necessary  condition  for  the  asymptotic  stability  of 
SLIP  for  an  arbitrary  leg  angle  trajectory  as  well  as  a 
sufficient  condition  for  its  instability.  These  conditions 
are  formulated  as  an  exact  algebraic  expression  despite 
the  non-integrability  of  the  SLIP  system.  Hence  leg 
recirculation  strategies  that  violate  the  above  condi¬ 
tion  can  be  discarded  without  recourse  to  cumbersome 
numerical  simulations.  We  also  use  the  closed  form  ex- 
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pressions  to  characterize  the  “cost”  of  sensing  required 
for  the  imposition  of  “fast”  transients  in  a  variety  of 
2  DOF  SLIP  models  that  have  appeared  in  the  recent 
legged  locomotion  literature. 

We  finally  apply  this  analysis  to  a  particular  3  DOF 
SLIP  model  with  pitching  dynamics  and  a  RHex-like 
leg  recirculation  strategy  that  satisfies  the  necessary 
condition  for  asymptotic  stability  in  certain  parame¬ 
ter  regions.  An  accompanying  numerical  study  shows 
that  this  model  captures  the  salient  aspects  of  the 
steady  state  dynamics  of  the  robot  RHex  [1]  (sim¬ 
ulated  in  SimSect  [7])  and  accurately  predicts  the 
robot’s  stability  properties. 

This  analysis  provides  for  the  first  time  a  partial  ex¬ 
planation  for  the  surprisingly  stable  behavior  observed 
empirically  in  the  robot  RHex.  It  also  paves  the 
way  for  a  more  principled  investigation  of  detailed, 
biologically-motivated  leg  placement  strategies  in  the 
LLS  model  [11]  which  captures  many  aspects  of  cock¬ 
roach  locomotion  [12]. 
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A  Time  reversal  symmetry  of 
RHex-like  leg  angle  trajecto¬ 
ries 

In  this  appendix  we  apply  the  condition  in  Lemma 
2  to  a  particular  family  of  leg  recirculation  schemes, 
thus  proving  the  involutive  nature  of  the  correspond¬ 
ing  time  reversed  flow  map.  In  particular,  we  prove 
that  the  solution  Itd{x o)  of  the  threshold  equation 
(19)  for  S2  =  G  o  f2  at  Jo  £  T  for  a  particular  leg 
angle  trajectory  also  solves  the  threshold  equation  at 
S%(x 0).  We  focus  on  the  family  of  leg  angle  trajecto¬ 
ries 

( j)(t ,  £0)  =  a(t)  +  k  (arccos(zo)  +  0o)  —  0o  —  0o 1  (58) 

where  a(t)  is  an  arbitrary  analytic  function  of  time. 
This  family  has  the  form  of  the  RHex-like  recirculation 
strategy  (55).  The  threshold  function  h2(fo(xo),xo,t) 
for  a  3  DOF  SLIP  model  reads 

h2(f2(x0),x0,t)  =  z(t)  -cos(<fi(t,x0))  ■  (59) 


Then  using  G(:ro)  : 

Gof*™(x0) 


(zo,-0 o,-z0,Oo)T  and 

^0  +  Z0tTD  - 
.  ~  (00  +  00 tTDJ 

—  (io  —  trn) 


V  00  J 

the  threshold  function  in  Lemma  2  reads 


(60) 


h2{G(x0),G  o  f*TD(x0),tTD)  = 
zo  —  cos  (a(tTi>)  ^  0o tm  +  fcarccos(zo  +  zotro— 

—  (00  +  OotTD){k  -  1))  =  0 

For  a  solution  of  this  equation  with  the  leg  recircu¬ 
lating  only  once  during  flight  <f>(tTD,x 0)  £  (|7r, 27t). 
This  must  be  taken  into  account  when  inverting  the 
cosine: 

arccos(zo)  =  -(a(tTD)  -  0otTD  + 

t2 

k  arccos(^o  +  zqItd - ^) 

—  (0o  +  0o  tTD){k  —  l)^j  +  2tt 

<=>  cos  (. k  arccos  (z^td)))  —  cos  (  arccos(z0)  +  a{txD ) 

—  00  tTD  —  (00  +  00 tTD){k  —  1)^  =  0 

O1  z{tTD)  -  COS  (<t>(tTD,X 0))  =  0 

with  Xo)  as  in  (58).  For  k  =  1  this  equation  is 
equal  to  the  original  threshold  equation  (59)  at  xo  and 
hence  tmix 0)  also  solves  the  threshold  equation  for 
S2(x 0).  Assuming  that  trn( Xo)  is  also  the  minimal 
solution  of  the  threshold  equation  at  S^^o)  f°r  ail 
xq  €  X,  we  can  conclude  that  the  time  reversed  flow 
map  of  the  flight  phase  with  a  leg  angle  trajectory 
defined  by  (58)  with  k  =  1  is  an  involution  on  Xh2  = 

X .  According  to  Lemma  3  this  means  that  F2oG  =  R2 
is  also  an  involution.  Then  \det(DxF2(x)\  =  1  at  a 
fixed  point  x. 


B  Equivalence  of  apex  and 
touchdown  Poincare  sections 

In  section  3.3.1  it  was  noted  that  in  [13,  15]  one¬ 
dimensional  Poincare  maps  characterized  by  the  apex 
event  during  flight  phase  were  used  to  illustrate 
the  asymptotic  behavior  of  the  constant  leg  touch¬ 
down,  leg  retraction,  and  “optimized  selfstabiliza¬ 
tion”  strategies  for  the  2  DOF  SLIP  model.  On  the 
other  hand  straightforward  counting  of  dimensions 
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shows  that  the  Poincare  section  of  a  2-dimensional 
SLIP  model  should  be  2-dimensional:  the  dimension¬ 
less  phase  space  X  :=  R  x  R+  x  R2  with  elements 
x  =  (y,z,y,z)T  is  four-dimensional;  conservation  of 
energy  E{x)  =  Eq  and  the  definition  of  the  Poincare 
section  V  :=  {x  €  X  :  p(x)  =  0,y  <  0}  should  re¬ 
duce  the  dimension  by  two.  For  the  Poincare  section 
denoting  the  touchdown  event,  p(x)  =  y/y2  +  z2  —  1 
whereas  for  the  Poincare  section  denoting  the  apex 
event,  pA{x)  =  z.  Using  conservation  of  energy  to 
eliminate  y,  the  reduced  Poincare  sections  can  then 
be  parametrized  by  x  =  (z,  z)T  €  X  for  the  touch¬ 
down  event  and  by  xa  =  (yA,  za)t  €  Xa  for  the  apex 
event.  While  for  some  singular  leg  placement  strate¬ 
gies  the  reduction  to  a  1-dimensional  Poincare  section 
at  the  touchdown  event  is  obvious,  e.  g.  for  the  con¬ 
stant  leg  touchdown  angle  strategy  illustrated  in  Fig. 
2  a),  there  is  a  priori  no  reason  why  the  apex  Poincare 
section  should  only  be  parametrized  by  one  variable. 
In  order  to  illustrate  the  equivalence  of  the  discrete  dy¬ 
namical  systems  defined  by  the  two  different  Poincare 
sections,  we  give  an  explicit  coordinate  transformation 
between  the  coordinates  of  the  two  reduced  Poincare 
sections  X  and  Xa-  Because  of  the  reset  of  the  y- 
coordinate  at  touchdown,  this  coordinate  transforma¬ 
tion  Ta  relates  the  touchdown  variables  x  to  the  next 
apex  variables  xa  and  not  to  the  previous  apex  vari¬ 
ables.  It  has  the  form 


Ta'-X  -  XA 


where  yA  =  \/l  -  z2LO  +  (^/2{E0  -  z)  -  z 2)  zLO  and 

za  =  ZLO  +  and  (zLO,  Zlo)  =  GoRi(x).  In  the 
notation  of  the  controlled  plant  model  (2),  the  control 
inputs  -  the  apex  to  touchdown  time  t-A  =  uA  and 
liftoff  to  touchdown  time  trn  =  u  -  are  related  by 

UA(k)  =  u(k)  -  ZLo(k)  .  (61) 

The  difference  between  the  two  parametrizations 
arises  in  the  structure  of  the  controlled  plant  model 
maps  A:  the  apex  controlled  plant  model  map  Aa 
decouples  into  two  separate  maps  Az  and  Ay  that  are 
independent  of  yA(k)  because  of  the  y  coordinate  reset 
at  touchdown:14 

(■ zA(k+l),yA(k  +  l))T  =  A((zA(k),yA(k))T  ,uA(k)) 
zA(k  + 1)  =  Az(zA(k),uA(k )) 
yA{k  + 1)  =  Ay(zA(k),uA(k)) 

14The  formal  expressions  for  Az  and  Ay  can  easily  be  derived 
and  are  not  given  here. 


Hence  the  only  way  that  yA(k)  can  enter  the  apex 
return  map  Ra  is  through  feedback:  uA(k)  =  tA(k)  = 
tA(zA(k),yA(k)).  Omitting  the  variable  yA(k)  which 
denotes  the  horizontal  distance  between  the  toe  pivot 
and  the  flight  phase  apex,  a  one-dimensional  return 
map  zA(k+  1)  =  Rz(zA(k))  =  Az(zA(k),tA(zA(k))) 
results.  For  the  3  DOF  SLIP  model  with  pitching,  the 
apex  return  map  without  yA  (k)  dependence  reduces 
the  dimension  from  4  to  3. 

This  explains  why  the  apex  Poincare  section  is  a  con¬ 
venient  parametrization  if  feedback  is  restricted  to  a 
subset  of  the  Poincare  section  coordinates.  However, 
the  touchdown  Poincare  section  seems  to  be  a  more 
natural  choice  for  the  description  of  physical  systems, 
since  the  touchdown  event  is  clearly  easier  to  detect 
than  the  apex  event,  which  requires  velocity  sensing. 

C  Invariant  tori  near  a  fixed 
point  of  2  DOF  SLIP  models 

In  this  section  we  establish  criteria  for  the  neutral  sta¬ 
bility  of  fixed  points  of  the  return  map  R  of  a  legged 
locomotion  model.  The  closed  circles  in  Fig.  3d)  sug¬ 
gest  the  existence  of  1-dimensional  f?-invariant  tori, 
on  which  R  acts  quasi-periodically.  This  is  reminis¬ 
cent  of  area-preserving  mappings  which  can  possess 
KAM-tori  (see  [35]  and  references  therein);  however, 
as  indicated  in  the  determinant  contour  plot  for  the 
stance  phase  alone  (Fig.  4),  area  is  in  general  not  pre¬ 
served  in  a  neighborhood  of  the  fixed  point  of  R ,  unless 
the  leg  placement  policy  is  designed  to  exactly  com¬ 
pensate  for  the  determinant  deviations  of  the  stance 
phase.  It  is  well  known,  on  the  other  hand,  that  re¬ 
versible  dynamical  systems  can  mimic  the  behavior  of 
Hamiltonian  systems  in  the  sense  that  they  can  also 
exhibit  KAM-tori  ([38],  for  a  review  see  [28]).  We  will 
show  how,  under  certain  assumptions,  a  theorem  on 
reversible  mappings  [36]  can  be  applied  to  establish 
the  existence  of  /t-invariant  tori  in  a  neighborhood  of 
a  fixed  point. 

C.l  A  theorem  on  invariant  tori  near 
a  fixed  point  of  reversible  diffeo- 
morphisms 

Before  stating  the  main  theorem,  several  definitions 
and  a  lemma  must  be  provided: 

Definition  4  (Involution  of  type  ( p,q ))  Let  x  £ 

Rw  be  a  fixed  point  of  the  involution:  G( x)  =  x.  An 
involution  is  said  to  be  of  type  (p,  q)  with  p  +  q  =  N 
at  x  if  the  characteristic  polynomial  of  the  Jacobian 
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of  G  at  x  reads  (— l)Ar(A  +  1)P(A  —  l)q .  This  is  the 
general  form  of  the  characteristic  polynomial  at  the 
fixed  point,  since  any  involution  can  be  written  in  a 
neighborhood  of  its  fixed  point  as  a  partial  reflection 
[39]. 

Definition  5  (Reversible  diffeomorphism)  A 

diffeomorphism  R  :  Rw  — >  Rw  is  called  reversible  with 
respect  to  the  involution  G  if  G  o  Ro  G  =  R~l . 

Lemma  4  (Composition  of  involutions)  The 

composition  R  =  R2  o  R\  of  two  involutions  R±  and 
f?2  is  reversible  with  respect  to  each  of  them,  i.e. 
R\  o  R  o  R1  o  R  =  id  =  f?2  0  R  o  R2  o  R-  Likewise,  a 
diffeomorphism  R  that  is  reversible  with  respect  to  the 
involution  R2  can  be  written  as  R  =  f?2  0  Ri  where 
R\  is  another  involution  [40]. 

Definition  6  (Symmetric  fixed  point  of  a  re¬ 
versible  diffeomorphism)  By  Lemma  4  a  diffeo¬ 
morphism  R  reversible  with  respect  to  the  involution 
i?2  can  be  written  as  R  =  f?2  0  Ri-  A  fixed  point 
x  £  Rw  of  R  is  called  symmetric  if  it  is  also  a  fixed 
point  of  R2  [28]. 

The  reduced  Poincare  map  for  SLIP  models  in  this 
paper  was  factorized  as  R  =  R2  °  R\  (24).  If  Ri  and 
i?2  are  involutions,  then  the  following  abridged  version 
of  Theorem  2.9  in  [36,  pages  147-152]  can  be  applied: 

Theorem  4  (Invariant  tori  near  a  fixed  point 
of  a  reversible  diffeomorphism  (M.B.  Sevryuk, 
1986))  Let  R  and  R\  be  diffeomorphisms  R,Ri  : 
R^  — »  M.N ,  analytic  in  a  neighborhood  of  a  common 
fixed  point  x  £  Iw  and  let  R  be  reversible  with  respect 
to  R\.  Assume  that  the  eigenvalues  {Ai,A;}i=1  ,...,at/2 
of  the  Jacobian  at  the  fixed  point  DxR(x)  satisfy  A ,  £ 
S1^— 1,1}  and  {Aj}j=1)..  ,,n/2  are  pairwise  distinct.  In 
addition  assume  that  R  is  nondegenerate,  i.  e.  31  £  N 
such  that  R  £  T}  (for  a  definition  of  T}  see  [36]). 
Then  the  following  holds: 

a)  In  any  neighborhood  of  x  £  Rw  there  exist  N  12- 
dimensional  tori  invariant  under  R  and  R±.  The 
action  of  R  on  these  tori  is  quasiperiodic,  and  the 
frequencies  of  this  action  are  constant  on  those 
tori. 

b)  There  exist  neighborhoods  Oe  of  x  £  Rw 

(lime^o  diam((De)  =  0,  Oei  C  Oe2  if  e  1  <  €2) 
such  that  lime^o  =  1  where  Qe  denotes 

IIlGS(v/g  I 

the  union  of  invariant  tori  in  Oe. 

c)  R\  is  an  involution  of  type  (N/2,  N/2). 


C.2  Application  to  2  DOF  SLIP  mod¬ 
els 

We  will  now  argue  that  this  theorem  can  be  applied 
to  the  2  DOF  SLIP  model  with  a  RHex-like  leg  recir¬ 
culation  (49)  with  k  =  1  as  suggested  by  Fig.  3d). 
The  recirculation  strategy  (49)  is  clearly  of  the  form 
(58),  hence  R2  is  an  involution  by  the  result  of  Ap¬ 
pendix  A.  In  Section3.1.1  it  was  shown  that  the  par¬ 
tial  stance  return  map  R\  is  also  an  involution.  Next 
we  need  to  show  that  Ri  are  analytic  at  the  fixed  point 
x  «  (0.8772, -0.0764)t: 

Analyticity  of  the  stance  phase  return  map  fac¬ 
tor  In  Section  2.4.2  the  analyticity  of  the  stance 
phase  flow  ft1<-x°')  Was  established.  The  correspond¬ 
ing  threshold  function  hi  (15)  is  analytic  in  x(t)  and 
Xo  if  C  7^  0-  By  the  implicit  function  theorem,  i  1  (3?o) 

will  be  analytic  as  long  as  A  (c  (7i(^o)))  It^xo)  ± 
0.  At  the  fixed  point  x,  (  =  1  7^  0  and 

|  (C  lt=ti(x)  =  yy+zz  «  0.6829  0.  Hence 

tffx)  is  analytic  at  x.  Since  the  composition  of  ana¬ 
lytic  functions  is  analytic,  Fi  =  fl'^fx)  and  S±  are 
analytic  at  x,  and  R\  is  also  analytic  at  x. 

Analyticity  of  the  flight  phase  return  map  fac¬ 
tor  In  Section  2.4.2  the  flight  phase  flow  f^x^  was 
seen  to  be  analytic.  We  focus  on  the  leg  angle 
trajectory  (49).  The  corresponding  threshold  func¬ 
tion  hi  (19)  is  analytic  in  xq  and  t  if  0  <  zo  < 
1.  By  the  implicit  function  theorem,  fg^’o)  wiU  be 
analytic  as  long  as  ^  (M/K^o),  %o,  t))  |t=tl(x0)  7^ 
0.  At  the  fixed  point  x,  0  <  z  «  0.8772  < 
1  and  TtJ.h2(f^G(x)),G(x),t))  |t=t2(G(x))  =  ~z  + 
sin(w(— 2z)  +  arccos (z)  +  tt)co  «  —6.6551  7^  0.  Hence 
t,2{G(x))  is  analytic  at  G(x).  Then  the  composition 
F2  =  f*2(-G(-x^  (G(x))  is  analytic  at  G(x ),  and  S2  and 
i?2  are  analytic  at  x. 

If  Ri  are  analytic,  then  the  composition  R  is  also  ana¬ 
lytic  in  x  and  by  Lemma  4  in  Appendix  C.l  the  return 
map  R  is  reversible  with  respect  to  both  R2  and  Ii\ . 
The  numerically  determined  eigenvalue  of  R  :  R2  — > 
R2  at  the  fixed  point  x  in  Fig.  3d  is  Ai  =  —0.6956  + 
i0.7185  £  S1  \  {—1,1}.  The  nondegeneracy  condition 
cannot  be  verified  rigorously  due  to  the  nonintegrable 
nature  of  the  return  map,  but  is  assumed  to  hold 
since  degenerate  diffeomorphisms  are  exceptional  in 
the  sense  that  they  constitute  a  variety  of  codimension 
one  [41].  Then  the  theorem  predicts  one-dimensional 
tori  around  the  fixed  point  which  are  invariant  under 
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R  and  R\  and  i?2,  two  of  which  are  plotted  in  Fig. 
3d.  The  quasiperiodicity  of  the  action  of  R  is  corrob¬ 
orated  by  the  numerically  determined  trajectory  of  R. 
Theorem  4  also  predicts  that  R\  and  R2  are  involu¬ 
tions  of  type  (1,1).  For  the  flight  phase  partial  return 
map  f?2  this  was  established  in  Section  3.1.2.  For  the 
stance  phase  partial  return  map,  the  eigenvalues  of  the 
Jacobian  of  R\  at  the  fixed  point  x  were  numerically 
determined  to  be  «  1,-1. 

D  SLIP  fitting  protocol 

D.l  SimSect  simulations 

1.  A  simulation  must  successfully  complete  10  sec¬ 
onds  of  simulation  time  without  crashing,  i.e.  by 
maintaining  upright  forward  locomotion. 

2.  A  simulation  must  be  stable  over  the  last  20 
strides.  A  stride  is  a  part  of  the  trajectory  be¬ 
tween  two  isolated  maxima  of  the  2-component 
of  the  COM-trajectory.  Stability  is  measured  in 
terms  of  the  maximal  deviation  of  averaged  linear 
and  angular  velocities  at  the  minimum  of  the  2- 
component  of  the  COM-trajectory  for  the  last  20 
individual  strides  with  respect  to  the  respective 
averaged  quantities  over  those  strides. 

3.  The  second  to  last  stride  of  the  simulation  is  se¬ 
lected  and  is  required  to  exhibit  a  flight  phase 
at  the  isolated  maxima,  i.  e.  none  of  SimSect’s 
six  legs  touch  the  ground  over  a  finite  amount  of 
time  around  these  maxima.  In  addition,  at  the 
the  minimum  of  the  2-component  of  the  COM- 
trajectory  all  three  legs  of  the  stance  tripod  are 
required  to  touch  the  ground,  whereas  all  three 
legs  of  the  flight  tripod  must  be  off  the  ground. 

Once  a  SimSect  simulation  has  passed  all  of  the 
above  conditions,  it  is  called  stable.  Then  a  3- 
DOF  SLIP  model  is  fit  to  the  simulation  data 
of  the  selected  stride.  The  simulation  data  re¬ 
quired  for  fitting  is  given  by  time  series  vec¬ 
tors  f/Simi  2/Sim,  2/Sim,  ^Sim,  ^Sim,  ^Sim,  ^Sim,  ^Sim,  ^Sim 
which  form  the  part  of  the  stride  trajectory  where  all 
three  legs  of  the  stance  tripod  of  the  SimSect  model 
are  on  the  ground  whereas  all  three  legs  of  the  flight 
tripod  are  in  the  air  (“stance  phase”).15  Here  the  y- 
coordinate  denotes  the  forward  position  of  SimSect’s 
COM,  2  denotes  the  vertical  position  of  the  robot’s 

15  The  stance  tripod  of  RHex  and  the  SimSect  model  is  formed 
by  those  three  legs  that  are  simultaneously  in  the  slow  phase  of 
the  RHex  clock  controller.  The  flight  tripod  is  formed  by  the 
other  three  legs  that  are  in  the  fast  phase. 


COM,  and  9  denotes  the  pitch  angle  of  the  robot  in 
the  sagittal  plane. 


D.2  Fitting  procedure 

The  equations  of  motion  for  the  SLIP’s  stance  phase 
(13)  and  the  equation  for  the  total  conserved  energy 
in  terms  of  dimensional  variables  can  be  written  as 


my 
m(z  +  g) 


-dyV(y,z,6) 

-d~zV{y,z,6) 


$  =  -d§V(y,z,0)  (62) 

^ 0 2  =  Eo-  V(y,z,0) -mgz-  y(y2  +  22) 


with  V(y,z,9)  =  f(C~Co ^(l+cggfi+cg^O^+cw^2), 
C  =  \/(y  +  A y)2  +  22,  and  tp  =  arctan  ^+.A^ .  We 
want  to  determine  the  a  priori  unknown  parameters 
Cf  =  (k,  Co,  egg,  eg c^,  Ay,  E0)  by  fitting  the  numer¬ 
ical  data  of  a  single  stance  phase  of  a  SimSect  sim¬ 
ulation  to  the  equations  (62).  The  first  five  compo¬ 
nents  of  Cf  are  parameters  that  determine  the  SLIP 
potential  V.  The  sixth  component  Ay  resets  the  y- 
coordinate  origin  and  hence  determines  the  y-position 
of  the  foothold  of  the  fitted  virtual  SLIP  with  respect 
to  the  stance  data.  The  fitting  parameters  also  include 
the  total  energy  E0,  because  in  SimSect  the  total  en¬ 
ergy  is  not  constant  due  to  damping,  frictional  losses 
and  hip  motor  torques. 


In  order  to  determine  the  fitting  parameters  c/  a  non¬ 
linear  fitting  procedure  (using  Matlab’s  ‘lsqcurvefit’) 
is  employed  that  computes 


min \\Ffit(cf,x)  -  Xfit\\l  (63) 

Cf 

~  ~  ~  ^  2  z  2 

where  X  =  (f/Sim,  ^Sim,  $Sim,  Vsim  A  ^Sim)> 

&  fit  (urf/Sim,  fh(2gim  A  (/)  , /^Sim,  (//2)0gjm) , 

and  Ffitifif ,  x)  is  the  expression  obtained  by  inserting 
x  into  the  right-hand  side  of  the  equations  (62).  Once 
a  solution  c/  has  been  found,  the  quality  of  the  fit 
must  be  quantified. 


D.3  Fitting  error  assessment 

Instead  of  using  the  residual  (63),  which  lacks  an 
intuitive  physical  interpretation  and  does  not  repre¬ 
sent  an  error  measure  in  phase  space,  we  compute 
fitting  errors  as  in  [4,  37].  The  assessment  of  the 
quality  of  the  fit  proceeds  in  two  steps.  First,  a 
SLIP  simulation  over  the  same  period  of  time  as  the 
data  trajectory  is  run  with  the  fitted  value  of  c/. 
The  initial  conditions  are  taken  to  be  the  positions 
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and  velocities  of  the  data  trajectory  at  the  mini¬ 
mum  of  isim-16  Second,  the  resulting  SLIP  trajec¬ 
tories  ysup,  ysup,  «SLiP)  -Zslip,  #slip,  #slip  are  com¬ 
pared  to  the  data  trajectories  by  L2  percent  errors: 


Air,  =  100 


H-^Sim  —  -XSLIP||2 


US: 


im  1 1 2 


(64) 


Here,  X  G  {y,  y ,  z,  z,  9,  9}  and  1 1  •  1 12  is  the  standard  2- 
norm.  In  an  effort  to  simplify  the  assessment  of  the  fit¬ 
ting  error,  the  quality  of  the  fit  is  reported  as  two  num¬ 
bers  —  the  average  L2  percent  error  of  the  cartesian 
coordinates  A yzL2  =  ( AyL2+AzL2  +  AyL2+AzL2)/4 
and  the  average  L2  percent  error  of  the  pitch  coordi¬ 
nates  A  9L2  =  (A  Gl2  +  A9l2)/2. 
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